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ABSTRACT 

We present the MOdel for the Rise of GAlaxies aNd Active nuclei (morgana), 
a new code for the formation and evolution of galaxies and AGNs. Starting from the 
merger trees of dark matter halos and a model for the evolution of substructure within 
the halos, the complex physics of baryons is modeled with a set of state-of-the-art 
models that describe the mass, metal and energy flows between the various components 
(baryonic halo, bulge, disc) and phases (cold and hot gas, stars) of a galaxy. These flows 
are then numerically integrated to produce predictions for the evolution of galaxies. 
The processes of shock-heating and cooling, star formation, feedback, galactic winds 
and super-winds, accretion onto BHs and AGN feedback are described by new models. 
In particular, the evolution of the halo gas explicitly follows the thermal and kinetic 
energies of the hot and cold phases, while star formation and feedback follow the results 
of the multi-phase model by Monaco (2004a) . The increased level of sophistication of 
these models allows to move from a phenomenological description of gas physics, based 
on simple scalings with the depth of the DM halo potential, toward a fully physically 
motivated one. We deem that this is fully justified by the level of maturity and rough 
convergence reached by the latest versions of numerical and semi-analytic models of 
galaxy formation. The comparison of the predictions of MORGANA with a basic set 
of galactic data reveals from the one hand an overall rough agreement, and from the 
other hand highlights a number of well- or less-known problems: (i) producing the 
cutoff of the luminosity function requires to force the quenching of the late cooling 
flows by AGN feedback, (ii) the normalization of the Tully-Fisher relation of local 
spirals cannot be recovered unless the dark matter halos are assumed to have a very 
low concentration, (iii) the mass function of HI gas is not easily fitted at small masses, 
unless a similarly low concentration is assumed, (iv) there is an excess of small elliptical 
galaxies at z = 0. These discrepancies, more than the points of agreement with data, 
give important clues on the missing ingredients of galaxy formation. 
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1 INTRODUCTION 

Several datasets spanning a large range of distances and lu- 
minosities in most of our past-light cone are now constrain- 
ing the background cosmology and the properties of primor- 
dial fluctuations with a remarkable precision. Temperature 
fluctuations of the CMB (Spergel et al. 2003, 2006), distant 
Supernovae (SN, Perlmutter et al, 1999; Knop et al., 2003), 
the large-scale structure traced by galaxies (2dF, Percival 
et al. 2002; SDSS, Eisenstein et al. 2005), the statistics of 
microlensing (Refregier 2003) , the abundance and clustering 
of galaxy clusters (Rosati, Borgani & Norman 2002) and the 



statistics of the Lymana forest transmission (Viel, Haehnelt 
& Springel 2004) are now giving a consistent picture of our 
Universe, well described by the ACDM model with param- 
eters (Cl ,ClA,^b,h,as) ~ (0.3,0.7,0.04,0.7,0.8) (with the 
last quantity ranging from 0.75 to 09). 

This "concordance" model provides the initial condi- 
tions for the evolution of perturbations, responsible for the 
gathering of dark matter (DM) into bound and relaxed ha- 
los. These DM halos host most of the astro-physical pro- 
cesses that rule the formation of stars and galaxies from 
the primordial gas. However, while the initial conditions are 
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specified and the basic physical laws are known, the forma- 
tion of galaxies is still an open problem, due to the high level 
of non-linearity of the processes involved and to the wide 
coupling of scales, from the km scale of imploding cores of 
SNe to the Mpc scale of galaxy super-winds. Galaxy forma- 
tion is thus a problem of complexity. 

Numerical simulations aimed at resolving galaxies (see, 
e.g., Pearce et al. 2001; Weinberg, Hernquist & Katz 2002; 
Steinmetz & Navarro 2002; Mathis et al. 2002; Lia, Porti- 
nari & Carraro 2002; Recchi et al. 2002; Toft et al. 2002; 
Springel & Hernquist 2003; Tornatore et al. 2003; Gover- 
nato et al. 2004) are still struggling both to tame the full 
complexity of the problem and to reach a sufficient mass and 
spatial resolution to resolve the "microphysics" of the injec- 
tion of energy by SNe and accreting BHs. As a matter of 
fact, most simulations of single galaxies are limited to mass 
resolutions of ~ 10 5 Mq and space resolutions of ~ 1 kpc, so 
the microphysical level is still treated as "sub-grid physics" 
and inserted in the codes with the aid of simple effective 
models. The problem is even more severe when trying to 
introduce accreting BHs within a galaxy (Di Matteo et al. 
2003; Kazantzidis et al. 2005), something that many authors 
regard as the missing ingredient of galaxy formation. Con- 
sidering then that most constraints on galaxies and AGNs 
are of a statistical nature, so that thousands if not hun- 
dreds of thousands of galaxies need to be produced for each 
model, it is clear that, at variance with what has happened 
with DM, a straighforward solution of the galaxy formation 
problem with N-body simulations is beyond hope for years 
to come. 

A quicker approach is given by the so-called "semi- 
analytic" galaxy formation models (White & Frenk 1991; 
Kauffmann et al. 1999; Somerville & Primack 1999; Cole et 
al. 2000; Wu, Fabian & Nulsen 2000; Granato et al. 2001; 
Hatton et al. 2003; Menci et al. 2004; Kang et al. 2005; Na- 
gashima et al. 2005; Cattaneo et al. 2006; De Lucia et al. 
2006; Bower et al. 2006), where all the processes that take 
place in the formation of galaxies are taken into account with 
simple approximated recipes. The main advantage of these 
models resides in the possibility of predicting the properties 
of whole galaxy populations in a short amount of computing 
time, thus making it possible to achieve a good sampling of 
the parameter space. However, it has been remarked that 
the (mostly phenomenological) recipes used in these mod- 
els have often a weak physical motivation and require the 
inclusion of free parameters that are difficult to constrain 
otherwise. As a consequence, the agreement between model 
and data, when achieved, may not shed much light on the 
process of galaxy formation. Moreover, given the intrinsic 
complexity of the problem, the models have often struggled 
to reproduce longly known pieces of evidence, like the shape 
of the luminosity function of galaxies, without convincingly 
showing their predictive power. This is highlighted by the 
difficulties of specific versions of semi-analytic models to re- 
produce some pieces of evidence, like the high-mass cutoff 
of the luminosity function (Benson et al. 2003), the nor- 
malization of the Tully-Fisher relation (Kauffmann et al. 
1999), the sub-mm counts (Baugh et al. 2005), the level 
of a-enhancement in ellipticals (Thomas et al. 2005; Na- 
gashima et al. 2005), the redshift distribution of K-band 
sources (Cimatti et al. 2002). From this point of view, it 



Figure 1. Outline of the main ingredients in the MORGANA code. 
Merger trees are taken from the PINOCCHIO tool (Monaco et al. 
2002a); galaxy mergers are treated following Taffoni et al. (2003); 
spectrophotometry is performed with GRASIL (Silva et al. 1998). 
The treatment of gas dynamics of baryons within the DM halos 
is the main argument of the present paper. 



is incorrect to claim that these models have too many free 
parameters, as it is clearly possible to falsify them. 

This paper is the first of a series devoted to presenting 
MORGANA, a new galaxy formation model which, in com- 
parison to the ones listed above, treats many of the physical 
processes at an increased level of sophistication. Figure [T] 
shows a general outline of this model, that could be applied 
to most similar models. The most relevant features of MOR- 
GANA are the following: 

• the evolution of the various components and phases of 
a galaxy is followed by integrating a differential system of 
equations along each branch of a merger tree, thus allowing 
for the most general (and non-linear) set of equations for 
mass, energy and metal flows; 

• the halo gas is treated as a multi-phase medium (hot 
gas, cooling flow, halo stars) and its evolution is described 
by a model that takes into account the thermal and kinetic 
energies of the hot and cold phases; this model treats cool- 
ing and infall as separated processes, takes into account the 
mass and energy injection from the galaxy to the halo (galac- 
tic winds) and from the halo to the Inter-Galactic Medium 
(IGM; galactic super- winds) ; 

• feedback and star formation are inserted following the 
results of the multi-phase model by Monaco (2004a, here- 
after M04; 2004b), plus an additional prescription for kinetic 
feedback; 

• accretion onto BHs and its feedback onto the galaxy 
(through radio jets and quasar-triggered galaxy winds) are 
built-in. 

The increased level of sophistication allows to move from a 
phenomenological description of gas physics, based on sim- 
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pie scalings with the depth of the DM halo potential, toward 
a fully physically motivated one. 

MORGANA has not been developed to construct a "the- 
ory of everything" for galaxies; we consider this model sim- 
ply as a precious tool (i) to understand the interplay and rel- 
ative importance of the various physical processes that take 
place in galaxy formation, (ii) to incorporate easily more 
processes that are thought to influence galaxies, so as to test 
their effects, (iii) to produce mock galaxy catalogues which 
reproduce particular selection criteria, in order to investi- 
gate the properties of galaxy surveys like sample variance 
and completeness level. The ultimate aim is to increase the 
predictive power of such galaxy formation models. 

On the other hand, the increase in the level of sophis- 
tication can only lead to an increase in the number of pa- 
rameters of the model. Many of these parameters are fixed 
either by observations (like the cosmological parameters) or 
purely gravitational N-body simulations (like the parame- 
ters related to galaxy mergings), while others can be fixed 
in principle with the aid of hydro simulations (like the pa- 
rameters related to gas cooling and super- winds). Besides, 
fundamental quantities like the energy of a single SN, the 
Initial Mass Function (IMF) of stars, the chemical yields or 
the parameters connected to the feedback and accretion onto 
BHs are very uncertain. This is again a problem of complex- 
ity, which will not be solved by neglecting important degrees 
of freedom. So, instead of decreasing the number of param- 
eters we aim to test our model versus a large number of ob- 
servational constraints, so as to fix all the parameters and 
propose predictions for upcoming and future observational 
campaigns. 

This paper presents the algorithm in full detail and 
some results that show the main potentialities, merits and 
problems of the model. The stellar mass function produced 
by this model has already been presented in Fontana et 
al. (2006). Other upcoming papers will address specific cos- 
mogonical topics, like the statistics of the AGN/quasar pop- 
ulation (Fontanot et al. 2006a), the assembly of the stellar 
mass of bright galaxies (Fontanot et al. 2006b) and the con- 
struction of the Stellar Diffuse Component in galaxy clus- 
ters (Monaco et al. 2006). The paper is organized as follows: 
Section [5] presents an outline of the model, describing its 
structure and the mass, energy and metal flows. The next 
sections, from [3] to 1111 describe in detail the treatment of 
all the physical processes inserted in the model; the reader 
willing to skip the details may go directly to section 1121 
Section [3] describes the merger trees of DM halos, section [4] 
the merging of galaxies, section [5] the physics of the halo 
gas phases, section [6] the formation of bulges and discs, sec- 
tionQthe modeling of stellar feedback, section [8]the produc- 
tion of metals, section [9] the accretion onto BHs. Section [lOl 
discusses the parameters involved in the model, while sec- 
tion [11] outlines the main post-processing phases. Section [121 
gives some basic results of the model and compares them 
to available observations, while section [131 gives the conclu- 
sions. Two appendices report details on the merging and 
destruction times of satellites and a very preliminary anal- 
ysis of stability with mass resolution. 

In this paper we use a "concordance" ACDM cosmology 
with parameters f2 = 0.3, Qa = 0.7, Q b = 0.044, cr 8 = 0.9, 
Ho = 70 km s _1 Mpc - . All physical quantities are scaled to 
this H value. The new WMAP results (Spergel et al. 2006), 



published when this paper was at a very advanced state of 
preparation, suggest erg ~ 0.75, and this could lead to some 
modest though significant changes in the predictions given 
in this paper. 



2 OUTLINE OF THE MODEL 

In this section we outline the model, describing the structure 
of the code and the system of equations for the mass, energy 
and metal flows that is integrated by the code. 

2.1 DM halos and galaxies 

The merger trees of DM halos are obtained using the PINOC- 
CHIO code (the details are given later in section [3. 1 P ■ This 
is not considered as an integrated part of the galaxy forma- 
tion code; any other code for generating merger trees can be 
used in its place. With respect to the Millennium Simulation 
(Springel et al. 2005), the PINOCCHIO trees do not give any 
information on the substructure of DM halos. 

Like N-body simulations, PINOCCHIO is based on realiza- 
tions of Gaussian initial conditions on a cosmological grid, 
so the mass resolution of merger trees is determined by the 
grid particle mass. Each DM halo that gets as massive as at 
least ten particle^] constitutes a starting branch of the tree; 
the time at which this happens is named appearance time. 
A galaxy is associated to each starting branch. When a DM 
halo merges with a larger one, it disappears as an individ- 
ual entity and becomes a substructure (or satellite) of the 
larger DM halo. The fate of substructures is followed using 
the model of Taffoni et al. (2003; section rO) : while the ex- 
ternal regions of the satellite DM halo are tidally stripped, 
its core, which contains its associated galaxy, survives for 
some time. A galaxy that is associated with an existent DM 
halo is named central; in general (though not always) the 
central galaxy is the largest in a DM halo. Galaxies asso- 
ciated with substructures are named satellites. Dynamical 
friction brings the orbit of the substructure toward the cen- 
tre, making it eventually merge with the main DM halo; 
when this happens the satellite galaxy merges with the cen- 
tral one. In some cases, the substructure and its associated 
galaxy can be destroyed by tides. Substructures can also 
merge between themselves, but this process is neglected in 
the present code. 

In general, DM halos will contain one central galaxy and 
several satellites, each associated with a substructure. When 
two DM halos merge, the central galaxy of the more massive 
DM halo will become the central galaxy of the merger, while 
the one of the smaller halo will become a satellite like the 
others. The bound between satellites belonging to the same 
substructure is assumed to be lost after the merger, i.e. we 

1 As shown in Monaco et al. (2002b), ten particles are not enough 
to reconstruct robustly a DM halo. On the other hand, in PINOC- 
CHIO the halos gather around the peaks of the inverse collapse 
time field, so the natural choice for the appearance time of a halo 
would be the collapse time of its first (peak) particle, correspond- 
ing to the creation of a 1-particle object. We deem that 10 parti- 
cles is a good compromise between the need for mass resolution 
and robustness. 
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do not allow for substructure of substructures. As a conse- 
quence, an existing substructure will always be associated 
with a single galaxy. 



2.2 Algorithm 

The flow chart of the algorithm is given in figure [2] The 
algorithm can be ideally divided into two main parts, repre- 
sented by the blocks on the left or right of the figure; the first 
part (in the left half of the flow chart) handles the merger 
trees and calls the numerical integrator for all the galaxies, 
the second part (in the right half) performs the integration. 
The algorithm proceeds as follows: 

• The merger trees are read from the PINOCCHIO output 
file. Then a loop on all the trees is started. 

• Each merger tree is subdivided into branches. A branch 
is defined as the evolution of a DM halo between two consec- 
utive mergings, be them major or minor; in other words, it 
corresponds to a time interval in which no new substructures 
are added to the DM halo. 

• The galaxies contained in the DM halo during the 
branch are looped on. The central galaxy is always addressed 
as the last one, so that the evolution of the hot halo gas as- 
sociated with it can take into account the energetic input of 
all the satellites. 

• For each satellite galaxy the dynamical friction, tidal 
stripping and tidal destruction times are computed (sec- 
tion gTj. 

• The time interval corresponding to the branch is sub- 
divided into smaller intervals, starting or ending at times 
isampicAt, where i samp i c is a sampling index and A t is usu- 
ally set to 0.1 Gyr. This is done to allow for a regular time 
sampling of baryonic variables. In general the branch ends 
are not integer multiples of At; for instance, if a galaxy is 
to be evolved from ii = 4.386 to ti = 4.728 Gyr, the branch 
will be subdivided into the intervals [4.386,4.4], [4.4,4.5], 
[4.5,4.6], [4.6,4.7] and [4.7,4.728]. 

• The evolution of the galaxy, the second main part of 
the code, is performed in the time interval by integrating 
a system of differential equations with a Runge-Kutta in- 
tegrator with adaptive time-steps (Press et al. 1992). The 
time-steps are chosen so as to have an accurate results to 
within 10~ 4 on the mass and energy flows described below. 
Of course, the integration of the galaxy stops at its merger 
or destruction time whenever these events occur. 

• The system of equations integrated by the code and the 
implemented physical processes are described below (sec- 
tion [2]4] a nd tables 1, 2 and 3). 

• The integration is interrupted whenever an instability 
is found; presently we take into account disc instabilities and 
quasars-triggered galaxy winds. 

• If relevant, the effect of such instabilities is applied at 
the end of the integration, and the loop on time intervals is 
closed. 

• At the end of the evolution of the galaxy in the branch, 
the predicted events of galaxy mergers, tidal stripping and 
tidal destruction, and the removal of the baryonic halo com- 
ponent of new satellites (described in section 12. 5[) are ap- 
plied. 

• The loops on galaxies and on branches are closed, and 



the results for all the galaxies in the tree are written on the 
output file. 

• The loop on trees is closed. 

2.3 Baryons 

The baryonic content of each galaxy is divided into three 
components, namely a halo, a bulge and a disc (figure [3}. 
Each component is made up by three phases, i.e. cold gas, 
hot gas and stars. The halo component of central galax- 
ies contains the virialized gas pervading the DM halo, cold 
gas associated to the cooling flow and halo stars. Satellite 
galaxies do not have any mass (or energy or metals) in their 
halo component, and this is justified by the fact that tidal 
stripping is very efficient in unbinding the halo component 
from satellites. (From the computational point of view, it 
is straightforward to relax this assumption, but we do not 
implement this feature here.) For the halo component, the 
code follows the mass and metal content of the three phases 
(cold gas, hot gas, stars), plus the thermal energy of the hot 
gas (that determines its temperature) and the kinetic energy 
of the cold gas (that determines its velocity dispersion). In 
total, 8 variables are associated to this component. 

Bulge and disc components are formally described in 
the code by the same variables. This approach was indeed 
undertaken as a follow-up of the M04 multi-phase model 
of a star-forming ISM, which is presented in some detail 
in section [7] However, a straightforward implementation of 
the M04 model, attempted in one of the first versions of 
the code, led to a number of unwelcome numerical compli- 
cations. We then decided to collect the main results of the 
M04 model and insert them as a set of simple recipes. This 
allows also to have a more transparent grasp of the effect 
of feedback in galaxy formation without losing the physical 
motivation of the multi-phase approach; this is a welcome 
feature in account of the lack of a widely accepted model 
for the evolution of a star-forming ISM. In summary, the 
bulge and disc components are described, in an effectively 
single-phase formalism, by a more limited set of 4 variables, 
namely the masses of stars and gas, and their metal masses. 

Together with the 16 variables (8+4+4) associated to 
the three components, 9 more variables are integrated by the 
code. The complete list of the 25 variables is the following: 

• mass, kinetic energy and metals of cold halo gas; 

• mass, thermal energy and metals of hot halo gas; 

• mass and metals of halo stars; 

• mass and metals of gas in bulge and disc; 

• mass and metals of stars in bulge and disc; 

• cooling radius of the hot halo gas; 

• mass, metals and kinetic energy of the cold gas ejected 
from DM halos as a super- wind; 

• mass, metals and thermal energy of the hot gas ejected 
from DM halos as a super- wind; 

• black hole mass; 

• black hole reservoir. 

The sampling of variables in the time grid defined above 
(section |2.2[) is performed as follows. At times i sa mpieAt the 
values of the 16 variables relative to the galaxy components 
are stored together with bulge and disc radii and velocities, 
black hole masses, punctual values for the star formation 
rates of bulge and disc and accretion rate onto the black 
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Figure 2. Flow chart of the MORGANA code; see text for all details. 



hole. Cooling radii, ejected matter and black-hole reservoirs 
are not sampled. With our choice of At = 0.1 Gyr the time 
sampling is adequate to describe old stars, but is too poor 
for young stars. For this reason we store the punctual val- 
ues of star formation in bulge and disc; in this way we can 
reconstruct, with a minimal amount of information, the re- 
cent star formation which gives rise to the massive stars. 
This procedure is presented and tested in Fontanot et al. 
(2006b) (see also section [TT]). 

2.4 Physical processes and flows 

The present paper is dedicated to a detailed description of 
all the physical processes included in the code. Here we give 
only a global view of these processes with references to the 
section of the paper where they are described. We also report 
the system of equations that is integrated by the code, and 
a list of the mass, energy and metal flows. 



Baryonic matter flows between the components and the 
phases as illustrated in figure [3] Within each component the 
three phases may exchange mass through the following flow 
terms: 

• evaporation of gas from the cold to the hot phase (M ev ); 

• cooling of gas from the hot to the cold phase (M co ); 

• star formation from the cold to the star phase (M s f); 

• restoration from the star to the hot phase (M ra ). 

As matter of fact, these flows are formally present in the 
code but they are never active at the same time. Indeed, 
in the halo component of central galaxies star formation 
(and then evaporation and restoration) are not active. In 
the bulge and disc components the hot phase is kept void by 
equating the restoration and cooling flow from the one hand, 
the evaporation and hot wind flow on the other hand; this 
is indicated in figure [3] by connecting, in the bulge and disc 
component, the arrows corresponding to the equated flows. 
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Figure 3. General scheme of mass flows in a model galaxy. The 
flows highlighted by dashed lines and gray colours are not used 
in the present version of the model. In particular, star formation 
(and then restoration and evaporation) in the halo is not active, 
while in bulge and disc components the evaporation and hot wind 
flows, as well as the restoration and cooling flows, are equated so 
as to leave the hot phase empty. 



In the following we will restrict ourselves to the flows that 
are used in the present formulation of the model, with the 
caveat that the other flows may be easily activated whenever 
required. 

Baryonic matter flows between the components as fol- 
lows: 

• primordial gas flows to the halo together with the ac- 
creted DM mass (cosmological infall, M cosm ); 

• cold gas infalls from the halo to the disc or bulge (M- m ); 

• cold and hot gas are expelled by the bulge or disc to 
the halo in a galactic wind (M cw and Mhw); 

• both hot and cold gas are allowed to leave the halo in 
a galactic super- wind (M hsw and Af csw ); 

• this expelled material is allowed to get back to the halo 
together with the cosmological infall; 

• winds from satellites are injected in the halo component 
of their central galaxy (M csat and iWhsat)- 

Mass conservation implies the following relations: 



M in ,H = M ift ,B + M in , D 
M CW , H = M CW , B + Af cw , D 
Mhw,H = Mhw.B + Mhw.D 



(1) 



As illustrated above (section 12. 2[) , in the cases of in- 
stabilities (disc instabilities and quasar-triggered galaxy 
winds), galaxy mergings, tidal stripping, tidal destruction 
and removal of halo baryonic component from new satel- 
lites, baryonic masses, energies and metals are moved be- 
tween components and galaxies outside the integration rou- 
tine. These flows are named in this paper external flows. It 
is worth noting that, as star formation in the halo is not 
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M B h = min(M vi sc , M B H:/*Ed) 

M r esv = Mi ow j — Mbh 



Cooling radius (equation 1281 1 



^cool = (M co ,h - M hw>ll )/(4wp g (r cool )r^ ool ) 



Table 1. The system of equations that is integrated by the code. 



allowed, stars flow to the halo component only through ex- 
ternal flows (namely tidal stripping and galaxy mergings). 

The system of equations integrated by the galaxy evo- 
lution code is reported in table 1. For all the equations, the 
suffixes c, h and s in the right hand sides denote the cold, 
hot and star phases. In all pedices the letters H, B, D and E 
following the comma denote the flows relative to the halo, 
bulge and disc components, and the matter ejected out of 
the DM halo by super-winds. 

Table 2 gives a list of all mass flows, with a quick expla- 
nation and a reference to the equation where they are de- 
fined; metal flows follow trivially mass flows (see section [8) 
with the exception of the newly generated metals, so for sake 
of brevity we only report these in the table. Finally, table 3 
gives the list of processes modeled in this paper, with refer- 
ence to the relative section and the list of the related mass 
and energy flows. 
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Physical process implemented in the integration 


Section 


Flows modeled in that section 


Galaxy winds from satellite galaxies 


[231 


M csat , -Kcsat, M hBBut , -Ehsat 


Hydrostatic equilibrium for the hot halo phase 


ED 


- 


Shock heating of cosmological infalling gas 


EH 


^cosm i -Ecosm 


Radiative cooling of hot halo gas 


[531 


Mco.H, S COi H, ^co,H 


Infall of cold gas on the galaxy 


[531 


Min,H, -ftTin.H, M in ,B, M in ,D 


Galaxy super-winds 


1531 


A^hsw, -E-hsw: -Eadi A^ csw , fSTcsw, A" a d 


Disc structure 


If; 11 
10- -LI 




Bulge structure 


[631 


- 


Star formation and feedback in discs 


[731 


Mrf,D, Afrs,D, Afhw.Di Mcw.D, E hWjH , -K"w,H 


Hj r j_ • 1 r 1 1 I'll 

Star formation and feedback in bulges 


17.31 


'"sf.B, M rs,B. M hw,B, M cw,B, «hw,Hi K w,H 


Metal enrichment 


IS 


hw,H' yi,B' yi,D 


Accretion onto BHs 


[931 


Advise i MlowJ 


AGN feedback 


[931 


-Ehw,H 








Disc instabilities 


[631 




Quasar-triggered galaxy winds 


[931 




Decay of satellite orbits by dynamical friction 


1431 




Tidal stripping at the periastron 


1431 




Tidal destruction of satellites 


1431 




Galaxy mergers 


1431 




Stripping of halo component of satellites 


1231 





Table 3. List of physical processes implemented in the code, with reference to the relevant section, and relative mass, energy and metal 
flows. 



2.5 Super-winds and satellite galaxies 

The matter ejected out of the DM halo by super- winds is col- 
lected into the variables denoted in table 1 by the ",E" pedix. 
After each integration interval these variables are stored in 
a vector, so as to be re-accreted at a later time together with 
the cosmological infall term; this is explained in section 15751 
In the present version of the model, satellite galaxies do 
not retain their baryonic halo component. However, satellite 
galaxies continually produce galaxy winds as long as star 
formation is active. The related flows are given to the halo 
component of the main DM halo as follows. First, the (cold 
and hot) winds and superwinds flows are equated: 



M csw 


= M CWj H 


-^CSW 


= J^cw,H 


Mhsw 


= Afhw.H 


-E^hsw 


= -Ehw.H 



Then, at the end of the integration over the time interval 
the ejected matter is added to a vector aimed to contain the 
contribution of all the satellites to the halo component of the 
central galaxy. We call these vectors M csa t , Mt sat and so on. 
The central galaxy is always evolved as the last one; when 
this happens the content of the satellite vectors is injected 
to the halo phase at a rate equal to the total mass or energy 
divided by the At time interval: 

M csa ± = Afcaat/At 

itcsat = max(A" csa t, Mcsat Ki 2 is p/2)/Ai 

M hsat = M hsat /Ai (3) 

-Ehsat = -Ehsat / At 



It is worth noticing that the kinetic energy of the cold wind 
is recomputed using the velocity dispersion of the DM halo 
as long as this velocity is higher than the kinetic velocity of 
the cold ejected gas. 

2.6 Initial conditions 

At the appearance time all DM halos are assumed to be as 
large as 10 particlefl All the bayons present in these pri- 
mordial DM halos are assumed to be in the hot halo phase, 
whose thermal energy, acquired by gravitational shocks, is 
computed with the model described in section 15.21 An issue 
with this setting is the quick start of cooling in these initial 
halos. This is in part a numerical artifact due to the lack of 
sampling of the tree; if this were resolved with a smaller par- 
ticle mass, the DM halo would then possibly contain some 
heating source at that time. To limit this resolution effect, it 
is assumed that the halos have just suffered a major merger 
at their appearance time, so that the onset of cooling is 
delayed by a few sound-crossing times (see section 15.21 for 
details). On the other hand, cooling in small halos is ham- 
pered by the ionizing background, and this is implemented 
(following Benson et al. 2002) by quenching cooling in all 
halos with circular velocities smaller than 50 km s _1 (sec- 
tion [OJ. 

Overcooling is connected to the very general problem of 
the stability of model predictions with respect to mass reso- 
lution. Appendix |B] presents a very simple convergence test 
of the model; the main conclusion is that model predictions 
do not converge and that the 50 km s _1 cutoff motivated by 

2 Some DM halos grow larger than 10 particles by a merger, but 
this detail is neglected. 
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Flow Comment Reference 



Mass flows 

A/ co h cooling flow eq. 1231 

A^in,H infall from halo eq. 1311 

Afcsw cold super- wind eq. 1421 

Af cw h cold wind to halo eq. [T] 

Mcsat cold wind from satellites eq. \3\ 

M C osm cosmological infall eq. 1191 

A^hsw hot super-wind eq. 1381 

A^hw.H h°t wind to halo eq. [T] 

A^hsat hot wind from satellites eq. \3\ 

A^in,B infall to bulge eos. 1331 1351 

M a { b star formation in bulge eq. 1631 

A^rs.B restoration in bulge eq. 1631 

A^hw B hot wind from bulge eqs. 1651 1731 1741 

A^cw,B cold wind from bulge eos. 1711 1751 

M in>D infall to disc eqs. [341 l36l 

M s f d star formation in disc eq. 1591 

A^rs.D restoration in disc eq. 1591 

A^hw D hot wind from disc eq. 1591 

A^cw,D cold wind from disc eq. 1591 

Kinetic energy of cold halo gas 

^co,H energy of cooling flow eq. 1251 

^in,H energy lost by infall eq. 1321 

Kcsw energy lost by cold super-wind eq. 1431 

^w,H energy acquired by winds eq. 1611 l67l 1721 

A"csat energy acquired by satellites eq. \3\ 

K&d adiabatic expansion eq. 1441 

Thermal energy of hot halo gas 

shock-heating of infalling IGM eq. r2T?l 

^co,H cooling of hot gas eq. 1241 

^haw energy lost by hot super- wind eq. 1391 

£hw,H hot wind to halo eqs. [60l [67l l83l 

^hsat hot wind from satellites eq. \3\ 

E a d adiabatic expansion eq. 1401 

Newly generated metals 

b new metals in bulge cold gas eq. 1781 

vi D new me tais in disc cold gas eq. 1781 

H new metals injected to the halo eq. 1771 

BH flows 

A^iowj loss of angular momentum eq. 1791 

Advise viscous accretion rate eq. 1801 



Table 2. Mass and energy flows in the system of equations (ta- 
ble 1). Metal flows are not reported. 



the reionization guarantees that convergence is forced at low 
redshift and high masses. This opens a very delicate topic, 
which has rarely been addressed in the literature (see Hat- 
ton et al. 2003); we leave a deeper analysis and discussion 
to further work, and limit ourselves to pointing out which 
of the results presented here are more sensitive to mass res- 
olution. In any case, as the lack of convergence regards star 
formation at very high redshift and faint galaxies, we can 
safely conclude that the building of bright galaxies does not 
depend strongly on mass resolution. 




Figure 4. Example of ordering of the PINOCCHIO merger trees. 
Circles with numbers denote DM progenitors, DM halos are de- 
noted by their linking lists. In all the mergings shown in this 
picture, the DM halo coming from the left is more massive than 
that coming from the right. 



3 DM HALOS 

3.1 DM merger trees 

As mentioned above, we use the pinocchio code to gen- 
erate the merger trees. The use of PINOCCHIO is motivated 
by its ability to give, even with modest computer resources, 
an adequate description of the hierarchical formation of DM 
halos, in excellent agreement with the results of N-body sim- 
ulations (see also Zhao et al. 2003; Li, Mo & van der Bosch 
2005); for instance, the mass function of DM halos is re- 
covered to within a 5-10 per cent accuracy (Monaco et al. 
2002a), while the mass function of progenitors of DM ha- 
los in a specified mass interval is recovered within a 10-20 
per cent error (Taffoni et al. 2002). This code uses a scheme 
based on Lagrangian Perturbation Theory (Moutarde et al. 
1991; Buchert & Ehlers 1993): starting from a Gaussian 
density contrast field sampled on a grid (very much simi- 
larly to the initial conditions of an N-body simulation) and 
smoothed over a set of smoothing radii, the earliest collapse 
time (the time at which the first orbit crossing takes place) 
is computed for each particle. Collapsed particles are then 
gathered into DM halos, where each halo is seeded by a peak 
of the inverse collapse time field. This procedure allows a de- 
tailed reconstruction of the DM halos, with known positions, 
velocities and angular momenta, and of their merger trees. 
The PINOCCHIO merger trees are equivalent to those given by 
N-body simulations, with a further advantage (shared with 
the less accurate Extended Press & Schechter approach; see 
Bond et al. 1992; Lacey & Cole 1993) of a very fine time sam- 
pling that allows to track the merging times without being 
restricted to a fixed grid in time (or scale factor). Other 
notable differences with respect to N-body trees is the im- 
possibility of PINOCCHIO DM halos to decrease in mass, a 
condition which is not strictly valid for the N-body simula- 
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tions. At variance with extended Press & Schechter merger 
trees, PINOCCHIO allows for multiple mergers of DM halos. 

The format of the PINOCCHIO outputs 
(in the updated 2.1 version, available at 
http : //adlibitum. oats . inaf . it/monaco/pinocchio/ ) 
are such that only the output at the final redshift is 
needed to reconstruct the merger trees of a realization, and 
from them compute the galaxy properties at all times. In 
PINOCCHIO, each DM halo retains its identity even when 
it disappears by merging with another larger halo. The 
output files contain, for each DM halo that has ever existed 
with at least 10 particles, the following information: (i) 
ID number, (ii) ID of the halo it belongs to at the final 
redshift, (iii) linking list of the halos, (iv) ID of the halo it 
has merged with, (v) mass of the halo at the merging time, 

(vi) mass of the halo it has merged with (before merging) 

(vii) merging redshift, (viii) appearance redshift. Halos that 
exist at the final redshift have field (i) and (ii) equal, and 
merging redshift equal to -1; field (iv) is also set to -1, while 
field (v) contains the mass of the halo at the final redshift 
and field (vi) is 0. 

Similarly to GALFORM (Cole et al. 2000), the linking 
list provided by the 2.1 version of PINOCCHIO is organized in 
such a way that merged halos are accessed so as to preserves 
chronological order. This is illustrated in the example of fig- 
ure [4] Each time a halo appears from a peak, its linking list 
points to itself. Suppose now that halo 1 merges with the 
smaller halo 2 at z merge =10, and that both halos have no 
substructure. Then the linking list is updated so that halo 1 
points to halo 2 and vice- versa. At z mergc =5.3 halo 1 (which 
contains halo 2 as a substructure) merges with halo 3, which 
has no substructure. Then the last halo of the chain, halo 2, 
is linked to halo 3, and this back to halo 1. Halos 5 and 6 
merge at z m erge=7.5, and their fate is similar to halos 1 and 
2. Halo 5 then merges with the larger halo 4 at z mcrgo =4.1; 
in this case halo 4 does not link directly to halo 5, which is 



put at the end of the chain, but to halo 6. At z„ 



=2.1 ha- 



los 1 and 4 merge. In this case, the last element of the chain 
of halo 1 (halo 3) is linked to the second of the chain of halo 
4 (halo 6). The final sequence is then 1-2-3-6-5-4. In more 
general terms, the two groups are linked in the following or- 
der: first the chain of the surviving DM halo, then the chain 
of the disappearing one, with the first element (the main 
halo before merging) put as last. It is clear that, starting 
from the second element of the final chain, two subsequent 
events are always accessed in chronological order. This does 
not imply a strict chronological order of all the mergings: 
in our example halo 3 (z mC rgc=5.3) is accessed before halo 6 
(zmer g e=7.5), but the two event are on independent branches 
of the tree. 



3.2 DM halo properties 

The physical properties of the DM halos (which are not 
predicted by PINOCCHIO) are computed at each integra- 
tion time-step (see figure [2} as follows. The density pro- 
file of the DM halo is assumed to follow Navarro, Frenk & 
White (1996; hereafter NFW), according to which a halo 
with a virial radius ru is characterized by a scale radius 
r s and a concentration c n f w = ru/r s - Defining the quantity 
S c = 200c„ fw (l + Cnfw)/ (3(1 + c nfw ) ln(l + c nfw ) - c nfw ), the 
NFW profile of a halo at redshift z is: 



Pdm(t) = pc(z) 



(4) 



where p c {z) is the critical density at the redshift z and x = 
r/r a . The virial radius ru of a halo of circular velocity Vh 
is computed assuming that its average density its 200 times 
the critical density (see, e.g., Mo, Mao & White 1998): 



(5) 



Vh = (10M H GH(z)) 1/3 
ru = Vn/WH(z) 

where Mh is the total halo mas fjfl and H(z) is the Hubble 
constant at z. 

The concentration c n f w is computed, as a function of 
Mh and z, following Eke, Navarro & Steinmetz (2001). 
Given the mass Mh, redshift z and concentration c n f w 
we compute the gravitational binding energy of the halo 
(Z7 H = J p$d 3 r) as: 



U H = 



/ii 



S c In (1 + c llfw ) 



1+Cn 



(6) 



From it and the virial theorem the velocity dispersion of the 
halo (such that the total kinetic energy of halo particles is 
= — Un/2) can be computed as: 



V diap = y/-Un/Mn 



(7) 



If A is the spin parameter of the DM halo, its specific 
angular momentum is computed as: 



■/ii 



-2- = GMhAV-O^h/Mh . 
IVlv 



(8) 



The DM halo parameters (Mh, c n f w , Th, Vh, Uu, Ju) are 
re-computed at each time-step along the integration. 

The spin parameter A of DM halos is in principle pro- 
vided by PINOCCHIO, which is able to predicted the angular 
momenta of halos. The predicted momenta show some (mod- 
est) degree of correlation with those of simulated halos (the 
spin directions tend to be loosely aligned within 60°), and 
their statistics is reproduced at the cost of adding free pa- 
rameters (Monaco et al. 2002b). While even a modest corre- 
lation with the N-body solution is an advantage with respect 
to drawing random numbers, the complications involved in 
reconstructing the spin history of halos weight more than 
any possible practical advantage. As a consequence, we pre- 
fer to randomly assign a A value to each DM halo, drawing 
it from the log-normal distribution: 



flog a (log A) a! log A = 



1 



2tt9 



exp 



(log A - log Ap) 2 
26» 2 



(9) 



where we use values Ao = 0.05 and 8 — 0.3 (Cole et al. 2000 
use the slightly lower value of 0.23 for 9; we follow Monaco, 
Salucci and Danese 2000 using 0.3, which is a better fit to 
many spin distributions available in the literature and cited 
in that paper). As explained in section [6~T1 a variation of A 
with time (for instance at major mergers) creates problems 
with the model of disc structure. To avoid these problems, 
the A parameters are held constant during the evolution of 



3 In the following we denote by Mh the total halo mass, including 
DM and baryons, and implicitly assume that the sub-dominant 
baryons do not influence the mass profile. This assumption will 
be relaxed in the computation of galaxy discs, section 16.11 
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each halo. This does not imply a constant specific angular 
momentum for the halos, as both Mh and Uh in equation [8] 
change with time. 

As for DM satellites, it is assumed that their prop- 
erties (mass, scale radius, concentration) remain constant 
after their merging. Satellites are however subject to tidal 
stripping, as explained in section H7T1 Tidal stripping is only 
taken into account in two cases: to strip baryonic mass from 
discs and bulges, in case of extreme stripping, and to com- 
pute the merging time for substructure after a merger. In all 
the other cases the unstripped mass of the satellite is used. 



4 GALAXY MERGERS, DESTRUCTION AND 
STRIPPING 

4.1 Dynamical evolution of satellites 

The computation of merging times for satellites follows the 
model of Taffoni et al. (2003) . In the simplest case, two DM 
halos without substructure merge, so that the smaller halo 
becomes a satellite of the larger one. The properties of the 
two halos at the merging time are computed as explained 
in section 13.21 and the subsequent evolution of the main 
halo is neglected (up to the next major merger). The orbital 
parameters are extracted at randonfl from suitable distri- 
butions, in particular the eccentricity of the orbit (defined 
as e = J/J c i r , where J is the initial angular momentum of 
the orbit and J c i r the angular momentum of a circular orbit 
with the same energy) is extracted from a Gaussian distri- 
bution with mean 0.7 and variance 0.2, while the energy of 
the orbit, parameterized with x c (defined as x c = r c /rn, 
where r c is the radius of a circular orbit with the same en- 
ergy and ru the halo radius) is taken to be 0.5 for all orbits. 
These numbers are suggested by an analysis of the orbit of 
satellites in a high-resolution N-body simulation (Ghigna, 
private communication), and are slightly different from the 
results of Tormen (1997), who found a lower value (0.5) for 
e and did not publish the distribution of x c (he gave the 
distribution of the radius of the first periastron, but this 
quantity is affected by dynamical friction in a subtle way); 
however, we have verified that choosing Tormen's value does 
not induce significant differences in the results. The choice 
of x c = 0.5 is equivalent to the implicit choice of Somerville 
& Primack (1999), who use a flat distribution for e, while 
Cole et al. (2000) extract a combination of the two parame- 
ters, e 7S Xc, from a log- normal distribution and Kauffmann 
et al. (1999) extract e again from a flat distribution and use 

Xc = 1. 

Galaxy mergings are due to the decay of the orbit of 
their host DM satellite by dynamical friction; in this scheme 
galaxy satellites can only merge with their central galaxy. 
Tidal shocks can lead to a complete disruption of satellites; 
in these rather unlikely events all the (dark and baryonic) 
matter of the satellite is dispersed into the halo of the central 
galaxy. Taffoni et al. (2003) give fitting formulae, accurate to 



4 We have verified that, as for the angular momentum of DM 
halos, PINOCCHIO is able to estimate the orbital parameters of 
merging halos, with roughly correct statistics. Again, we prefer 
to extract these parameters from suitable distributions. 



the 15% level, for the merging and destruction times for sub- 
structures that take into account dynamical friction, mass 
loss by tidal stripping and tidal shocks. We use a slightly 
different version of those fitting formule, with updated pa- 
rameters that slightly improve the fits to simulations; these 
are given in appendix |A1 

Though the merger and disruption times of Taffoni et 
al. (2003) include a rather sophisticated treatment of tidal 
stripping, we implement this process at a very simple level. 
Tidal stripping is applied at the first periastron of the satel- 
lite. The tidal radius is computed as the radius at which the 
density of the unperturbed satellite is equal to the density 
of the main DM halo at the periastron. All the mass exter- 
nal to the tidal radius is then considered as unbound from 
the satellite. However, while the stripping radius is recorded, 
the halo mass is not updated but kept fixed in the following 
evolution. Indeed, disc structure (see section 16. 1[) is com- 
puted using the DM halo profile at of the satellite at its 
merger time (i.e. ignoring stripping), while the recomputa- 
tion of merging times after a major merger is done using the 
stripped mass. 

Tidal stripping at periastron affects also the fraction 
of stellar and gas mass of disc and bulge that lies beyond 
the tidal radius. Notice that, for simplicity, neither the DM 
halo nor the disc and bulge are assumed to be perturbed 
by this process beyond the decrease of mass. Within the 
model structure, stripping is applied outside the numerical 
integration, so it is an external flow (section |2.2|I . 

In the more general case of two substructured halos that 
merge, we distinguish between major and minor mergers as 
follows: 

Major merger of DM halos : M sat /M main > / hm m (10) 

(we recall that the main DM halo includes the satellite). 
N-body simulations show that when this condition, with 
/hmm — 0.2, is satisfied the perturbation induced by the 
satellite leads to a reshuffling of all the orbits. At a ma- 
jor merger we then re-extract the orbital parameters and 
re-compute the merger and destruction times for all the 
satellites of the main DM halo, as if they just entered the 
halo. As a consequence, some galaxies near to their merging 
time can be moved to a different orbit that does not lead to 
merging, while some other galaxies can suffer tidal stripping 
more than once. Clearly the re-computation of the merging 
and destruction times for a substructure may not be very 
accurate, especially for satellites that have suffered strong 
mass loss. In this case we keep the scale radius and concen- 
tration of the satellite fixed, but use (as mentioned above) 
the stripped mass to compute the merging and destruction 
times. As these times are rather long for small satellites, the 
final results is that the galaxies just don't merge and the 
accuracy of the prediction is not an issue. 

Minor mergers do not influence the evolution of the 
satellites of the main DM halo, but do affect the satellites 
of the smaller DM halo (going itself to become a satellite), 
for which there is no difference between minor and major 
merger. 

4.2 Galaxy merger trees 

The galaxy merger trees are constructed, analogously to the 
DM halo merger trees, by specifying for each galaxy (i) the 
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galaxy where its stars lie at the the final time, (ii) the merg- 
ing redshift, (iii) the galaxy it has merged with, (iv) a linking 
list for the merged galaxies, (v) and (vi) the masses of each 
pair of merged galaxies. Destroyed galaxies are recorded by 
assigning a negative value to field (i). The construction of 
galaxy trees is performed as follows: at each halo merger the 
merging and destruction times for the galaxies are computed 
(section 14. ip , then the galaxies are merged or destroyed at 
that time if the DM halo they belong to has not been in- 
volved in a major merger nor has become a satellite in the 
meantime. While multiple mergers are allowed by PINOC- 
CHIO, galaxy mergers are all binary. 

4.3 Galaxy mergers 

When two galaxies merge, their fate depends again on the 
ratio of their masses. Major mergers of galaxies are defined 
as: 



Major merger of galaxies : M sa ,t/M ccn > f g 



(11) 



where the parameter / gm m is suggested by simulations to 
take a value of 0.3 (Kauffmann et al. 1999; Cole et al. 2000). 
In this specific case baryonic masses are used (i.e. the mass in 
hot, cold and star phases of the bulge and halo components), 
and the central galaxy does not include the satellite, so this 
condition is similar to that of equation \W\ with a value of 
0.25. While the condition on DM halos can be computed 
directly from the PINOCCHIO trees and without running the 
galaxy evolution code, the condition of equation [11] must be 
computed at the merger time, after the evolution code has 
determined the baryonic galaxy masses. 

At minor galaxy mergers the whole satellite is added 
to the bulge, while the disc remains unaffected. This is at 
variance with Cole et al. (2000), that give the stars to the 
bulge and the gas to the disc. Their choice is however ques- 
tionable, as the dissipative matter is more likely to go to 
the bulge; this is why we prefer to give everything to the 
bulge. A more accurate solution of this problem will clearly 
depend on the orbit of the merger, but an implementation 
of these second-order effects would require a large set of 
N-body simulations to quantify and parameterize these de- 
pendences. In any case, our tests have revealed no strong 
difference between the two cases (see section 16.21 for more 
comments) . 

At a major galaxy merger all the gas and stars of the 
two merging galaxies are given to the bulge of the central 
one. Due to the shorter time-scale of star formation in bulges 
(see section [7} , a starburst is stimulated by the presence of 
gas in the bulge component. 

Stripping and galaxy destruction (which is an extreme 
stripping event) bring stars to the halo of the central galaxy. 
As star formation in the halo is inactive, this is the main way 
to have stars to the halo. We anticipate that this mechanism 
is not very effective in galaxy clusters, where only a few per 
cent of the stellar mass is stripped to the halo, at variance 
with the 10-40 per cent found in observed clusters (Feld- 
meier et al. 2003; Arnaboldi et al. 2003; Gal- Yam et al. 2003; 
Zibetti et al. 2005). Murante et al. (2004; 2006) have per- 
formed hydro simulations to address this problem, coming 
to the conclusion that the high fraction of cluster stars can 
be reproduced, but the main mechanism is not tidal strip- 
ping but violent relaxation in major mergers. Very recently, 



Monaco et al. (2006), based on MORGANA and on the N-body 
results of Murante et al. (2006), proposed that the construc- 
tion of the stellar diffuse component of galaxy clusters, i.e. 
the halo star phase, is related to the apparent lack of evolu- 
tion of the most massive galaxies since z = 1. They showed 
that the evolution driven by galaxy mergers alone is larger 
than what is observed in large galaxy samples, and that this 
discrepancy can be solved if a fair fraction of stars is scat- 
tered to the halo star component at each merger. However, 
to reproduce the observed increase of the stellar diffuse com- 
ponent with halo mass, the fraction of scattered stars must 
depend strongly on the properties of the DM halo and of 
the merging galaxies. This dependence should be provided 
by simulations, but this is a work in progress. In this paper 
we simply assume that a fraction /scatter of the star mass of 
the merging galaxies is scattered to the halo at each major 
merging. With this simple rule we are able to produce a low 
fraction of halo stars in galactic halos and a higher one in 
galaxy cluster. 

Interactions between satellites, like binary mergers, fly- 
bies (that stimulate star formation) or galaxy harassment 
(Moore et al. 1996) are not included at the moment. We 
know that these events can have an impact on the evolution 
of galaxies (Somerville & Primack 1999; Balland, Devriendt 
& Silk 2003; Menci et al. 2004) and AGNs (Cavaliere & Vit- 
torini 2000). We plan to introduce such events in the future. 



5 HALO GAS 

5.1 Equilibrium model for the hot phase 

The hot halo phase is assumed to be (i) spherical, (ii) in 
hydrostatic equilibrium in an NFW halo (Suto, Sasaki & 
Makino 1998), (iii) filling the volume from a cooling radius 
''cool to the virial radius th, (iv) pressure-balanced both at 
r CO oi and ru, (v) described by a polytropic equation of state 
with index 7 P , a parameter that is observed to take a value 
of about 1.2 (see, e.g., De Grandi & Molendi 2002); we use 
the best-fit value of 1.15. The equilibrium configuration of 
the hot halo gas is computed at each time step as described 
below, under the assumption that in the absence of major 
mergers the gas re-adjusts quasi-statically to the new equi- 
librium configuration. We do not assume that the thermal 
energy Eh of the hot halo phase is a fixed fraction of the 
virial energy; on the contrary, we follow its evolution through 
the equation given in table 1. This is done to incorporate the 
response of the hot halo phase to the thermal feedback from 
galaxies. Under the assumptions given above, the constraints 
on the profile are then the mass and thermal energy of the 
gas. 

The equation of hydrostatic equilibrium: 



dP g = c P 9 Mu(< r) 
dr r 2 



(12) 



(where P is the gas pressure, r is the radius, Mh(< r) is 
the halo masijfl within r and p g the gas density) is easily 



5 Gravity is supposed to be dominated by DM, so we use for 
Mh(< r) the NFW mass profile, as if the baryons were distributed 
as the DM. The error induced in this assumption is not likely to 
be relevant; a more sophisticated treatment would not allow to 
obtain an analytic solution for the gas profile. 
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solved in the case of an NFW density profile (equation 
The solution is: 



Pair) 
PAr) 



Pgo 1 - a 1 



= P g0 I 1 - a I 1 - 



ln(l + c nfw a:) 

Cnfw 

ln(l + c nfvr x) 

Cnfw 



l/(7p-l) 



7p/(7p-l) 



(13) 



= T g0 



( ln(l + c n f w x) 
f — a I 1 

Cnfw 



Defining the virial temperature of the halo as T v i r = 
HhotrripVc /3k (where /ihot = 4/(8 — 5Y — 6Z) is the mean 
molecular weight of the hot gas; Z = Mjf H / 'M^h is the 
metallicity of the hot halo gas and Y assumes the solar value 
of 0.25), r) = T g (r)/T v i T and rjo its extrapolation at r — 0, 
the constant a is defined as: 



7p 



1 3 



Cnfw ( 1 H~~ C n f w ) 



7p VO (1 + C nfw ) ln(f + C n fw) — C n fw 



(14) 



The constants Tgo, PgO and P s o are defined as the extrap- 
olation of the density and temperature profiles to r = 0, 
even though the gas is assumed to be present only beyond 
fcooi- As mentioned above, they are fixed by requiring that 
the total mass and energy of the gas correspond to M n ,H 
and -Eh.H- The first condition can be solved explicitly if the 
energy is specified. Calling I the integral: 



1(a) = 



l/r. 



l-O 1 - 



ln(l+t) 



(15) 



(where for simplicity we declare only the dependence on the 
a exponent) we have: 



PgO = 



M h 



4-Kri 



2W(7p - 1)) 



The second condition is 

6irkT g0 pg r 3 



x Z(7p/(7p - 1)) 



(16) 



(17) 



MhotWlp 

This equation cannot be solved explicitly, as the coefficient 
a depends on the energy itself through rjo. To find a solu- 
tion to these two equations it is necessary to use an iter- 
ative algorithm. As a consequence, the computation of the 
two integrals contained in equations [16] and [17] is the most 
time-consuming computation of the whole code. A dramatic 
speed-up (at the cost of a negligible error) is obtained by 
computing the integrals in a grid of values of a, r and 7 P ; 
the solution is then found by linearly interpolating the table. 

The function 1 — a(l— ln(l + c n f w a;)/cnfw) in eauationsfTBI 
becomes negative at large radii. In this case density, pres- 
sure and temperature are not defined. Usually this happens 
beyond the virial radius th, unless the central temperature 
is lower than the following limit: 



rio < 3 



7p 



1 



Cnfw - ln(l + C nf 



7p ln(l + Cnfw) — C n f„/(1 + Cnfw) 



(18) 



This condition can be met at high redshift, when c n f w -values 
are high. In this case no gas is assumed to be present beyond 
the point of zero density, so that the gas is bound to the inner 
part of the halo and its pressure at the virial radius is null. 

This model for the hot halo gas is very similar to that 
proposed by Ostriker, Bode & Babul (2005) to model the 



hot component of galaxy clusters, with two remarkable dif- 
ferences: first, they do not follow cooling, so the hot gas is 
present since r = 0, instead of r = r CO ol- Clearly this does 
not make any difference in cases like galaxy clusters without 
cool cores, where r coo i ~ 0. Second, they assume an external 
pressure, computed on the basis of a fiducial infall velocity 
of cold gas, and extrapolate the gas profile until its thermal 
pressure equals the external one. This implies that the hot 
gas can extend beyond the virial radius. Our choice is to 
remove this outlier gas, and this allows to describe galaxy 
super- winds (section |5.5[) . Clearly the two criteria should be 
equivalent in the case of a small cooling radius and a thermal 
energy roughly similar to the virial one; we plan to deepen 
this point and to compare the predicted properties of the 
hot halo gas with simulations in the future. 



5.2 Shock heating 

The equilibrium model does not specify the amount of ther- 
mal energy of the hot gas. This is acquired by the infalling 
gas through shocks. The cosmological infall mass flow is 
computed by linearly interpolating the DM halo mass be- 
tween the branch ends (whose distance in time is At) and 
assuming that a fraction Qt/Qo of that mass is in IGM: 



Oj, A Mb 
Th At 



(19) 



We then assume that this gas acquires an energy equal to 
/shock times that suggested by the virial theorem: 



E c 



= /shock (-0.5 Un)M aosm /Mn 



(20) 



where the binding energy of the halo (7h is given by equa- 
tion [6] The parameter /shock is suggested by hydro simu- 
lations to be slightly higher than 1 (Wu, Fabian & Nulsen 
2000). We adopt a value of 1.2. 

A similar heating is applied in the following external 
flows: 

• the hot gas contained in the DM halos at the appear- 
ance time (see section [2~6]l : 

• the hot halo gas of satellite DM halos at their merging 
time' 

• the hot halo gas at major mergers. 



In all these cases: 



Ei 



/shock(-0.5(7 H )Mh,H/MH 



(21) 



Hydro simulations suggest that any cold gas present 
in the halo is re-heated by shocks during a major merger. 
Accordingly, we allow shock-heating to affect also the cold 
halo phase at major mergers. This option is active in all the 
results presented here, but can be switched off on request. 

Cases (ii) and (iii) refer to gravitational heating due to 
merging events. This heating is of course not instantaneous; 
the energy is redistributed to the whole halo gas in a few 
crossing times. This behaviour is implemented by quenching 
cooling for a number n qucnc h of crossing times fh/Vh; after 
the quenching (which ends at some time t q ), the cooling flow 
is allowed to start gradually as exp{ — [rn/Vn(t — t q )] 2 }. The 
parameter n quenc h is very important to control the cooling 
of gas, especially at high redshift. The same quenching is 
applied to DM halos at their appearance time (section 12.61) : 



© RAS, MNRAS 000,[TJj35] 



MORGANA 13 



it amounts to assuming that the appearing halos have just 
formed by a major merger. 

Using hydro simulations, some authors (Kravtsov 2003; 
Keres et al. 2005) have recently pointed out that cold gas 
can flow directly to the core of small-mass halos at high 
redshift by radiating very quickly the energy acquired by 
shocks. Following this idea, Dekel & Birnboim (2006; see also 
Cattaneo et al. 2006) have proposed that this lack of shock 
heating is a likely responsible for the observed transition in 
the behaviour from dwarf to bright galaxies. Besides, Croton 
et al. (2006) implement this idea by equating the infall mass 
flow with the cosmological infall one in the infall-dominated 
halos. A similar view is taken by Bower et al. (2006), who 
suppose that a short cooling time prevents the formation of 
a hydrostatic hot halo phase. These considerations question 
the validity of the shock heating and hydrostatic equilibrium 
hypotheses in infall-dominated halos; however, these are also 
the cases where the infall of gas on the galaxy does not 
depend much on the cooling time, so we choose to retain 
these hypotheses as they are accurate in the cases where 
they are most relevant. 



5.3 Cooling 

The cooling radius is defined as the radius within which the 
hot halo gas has cooled down. In most semi-analytic models 
the hot gas profile is computed at a major merger; the time- 
dependent cooling radius is then computed as the radius at 
which the cooling time of a gas shell is equal to the time 
since the merger. In the present model the cooling radius is 
instead treated as a dynamical variable. This allows to re- 
compute the gas profile at each time-step, and to incorporate 
the heating effect of the hot wind coming from the central 
galaxy. 

The cooling rate of a shell of gas of width Ar at a radius 
r is computed as: 



AM co , H (r) 



AMi,, H (r) 
tcooi(r) 



(22) 



where AMh,H = 4irr 2 p g (r)Ar is the shell mass and 
tcooiM = 3kT g (r)p hot m p /2p g (r)A cool (T g (r)) is the cooling 
time at radius r. For the cooling function A coo i, we use the 
met allicity- dependent function tabulated by Sutherland & 
Dopita (1993). The cooling time depends on density and 
temperature, but the density dependence is by far stronger, 
both intrinsically and because the temperature profile is 
much shallower than the gradient profile. So the integra- 
tion in r can be performed by assuming T 9 (r) ~ T 9 (r coo i). 
The resulting cooling flow is: 



AL 



co,H 



47rr|p 9 o 



x T(2/( 7 , - 1)) 



(23) 



where X is defined in equation 1 151 In this equation the cool- 
ing time £ C ool,o is computed using p g o for the density (the 
dependence of density on radius is taken into account by the 
integrand), and T 9 (r coo i) for the temperature, as explained 
above. Analogously to the integrals of equations 1161 and 1171 
the integral in equation [23] is computed on a grid of param- 
eter values and then estimated by linear interpolation on 
the table. The rate of energy loss by cooling is computed 
analogously: 



h = 3fcT 9 (r cool )4^p 9 o x I(2/(7p _ 1)} 

2/ihot'Tlp £cool,0 

The cooled gas carries with it a kinetic energy: 
1 • 2 



(24) 



(25) 



where Vdi sp is the velocity dispersion of the DM halo defined 
in equation [7] 

When a heating source is present, these two terms be- 
have differently. While the energy radiated away by the hot 
gas at a given density and temperature does not change, 
the amount of cooled mass depends on how much of this 
energy is replaced by the heating source. We then compute 
the cooling time as: 



icool.O 



3kT g (r cool )p hot m p 



(26) 



2p 9 o(A coo l — Theat) 

This cooling time is used in equation[23]to compute the mass 
cooling flow. A negative value implies a net heating of the 
source, in which case the mass cooling flow M co ,u (but not 
Sco.h) is set to zero. 

The source of heating is the central galaxy, which ham- 
pers cooling through the hot wind energy flow, i?hw,H- This 
flow carries the energy produced by SNe, both in the bulge 
and in the disc, and by the AGN. Satellites instead are as- 
sumed to orbit on average in the external regions, so that 
the energy contributed by their winds is injected beyond 
the cooling radius and does not interact directly with the 
cooling flow. To compute the heating term it is necessary 
to specify how this heating is distributed. For simplicity we 
assume that heating and cooling affect the same gas mass 
47rraP 9 oI(2/(7 P — 1)), i.e. the inner shell at r coo i that is ef- 
fectively cooling. The heating term is then computed as: 



The 



-Bhw,H 



47rrIX(2/( 7p - 1)) ^ p g0 



(27) 



Once the cooling and heating sources are fixed, the evolution 
of the cooling radius is computed by inverting the usual 
relation, M co ,h = — 47rp 9 r^ ool dr coo i/dt, taking into account 
that the hot wind mass flow Mh w ,H is adding to the hot halo 
phase at the cooling radius: 



f"cool 



co,H 



4np g (r cool )r^ 



(28) 



In this way, the cooling radius decreases if the hot wind term 
overtakes the cooling term. 

Equation [28] shows that the cooling radius should not 
vanish. Moreover, in the case of very strong cooling flows the 
Runge-Kutta integrator may try some calls with r coo i > m, 
giving rise to numerical problems. The cooling radius is then 
forced to lie between a small value (taken to be 0.01 times 
the scale radius r s ) and 90 per cent of the virial radius th- 
This is done by gradually setting r CO oi to zero when the 
limits are approached. The presence of a small lower limit 
for r coo i does not influence the results, because of the flat 
density profile in the central regions. The upper limit can 
instead influence significantly the behaviour of cooling, but 
this happens when most of the hot gas has already cooled, a 
situation in which a precise modeling is not very important 
and the validity of the model itself is more doubtful (see the 
discussion in section [5~TT l . 

The existence of a "cooling hole" at the centre of the 



© RAS, MNRAS 000,[TH35] 



14 Monaco, Fontanot & Taffoni 



halo requires that, as we assumed in section 15.11 pressure 
is balanced at r CO ol- However, this assumption is clearly 
unphysical, as the cooled gas cannot give sufficient pres- 
sure support, so this assumption will work only as long as 
^cooi > c s , where c s is the sound speed of the hot gas. In 
general, the gas at r coo i will be pushed toward the centre by 
its pressure. This can be modeled very simply as follows: 



^cool — ^cool 



(29) 



We have noticed that this feature induces an excessive and 
unphysical degree of cooling in the later stages of evolu- 
tion, when, due to the decreasing temperature gradient, the 
residual thermal energy of the gas gets lower than the virial 
energy. To avoid it we consider the sound speed term only 
when the gas thermal energy i5h,H is higher than the virial 
value. This term is included on request, and is used in the 
results presented here. 

Finally, after reionization the ionizing background is 
likely to prevent the cooling of any halo whose circular ve- 
locity is smaller than ~50 km s _1 . As suggested by Benson 
et al. (2002), to mimic the effect of the UV background we 
simply suppress any cooling in all halos with V c < 50 km 



The corresponding loss of kinetic energy is: 



^dyn^dyn (^"cool ) 



(32) 



The infalling cold gas is distributed between the disc 
and the bulge as follows. As a first option, all the infalling 
gas is given to the disc, under the assumption that it has 
the same specific angular momentum as the DM halo: 



M in , B 
M in , D 





M in , H 



(33) 
(34) 



In presence of a significant or dominant bulge the formation 
of such a disc by infall implies that the bulge has no influ- 
ence on it, even when a large fraction of it is embedded in 
the bulge. This is a rather strong assumption, as the hot 
pressurized phase pervading the bulge fsection l7.3[l can lead 
to significant loss of angular momentum of the gas by ram 
pressure. Then, as a second option we let gas infall on the 
bulge by a fraction equal to the mass of the disc contained 
within the half-mass radius of the bulge: 



5.4 Infall 

The dynamical time of the halo at a radius r is defined as 
the time required by a mass particle to free-fall to the centre: 



tdyxx(r) 



( 3<5 C 



y/2V c V 200c n 



-1/2 



(30) 



r I r s 



ln(l + y) ln(l + r/r s ) 



r/r 



'1/2 



dy 



The cold phase is unstable to collapse and flows to the cen- 
tral galaxy. Starting with White & Frenk (1991), many semi- 
analytic codes, at variance with MORGANA, unify the pro- 
cesses of cooling and infall by computing an infall radius for 
the gas as the radius at which the infall time of a gas shell 
is equal to the time since last major merger, then using the 
smallest between the cooling and infall radii to compute the 
cooling flow. This choice implies an assumption of no dif- 
ference between the hot halo gas and the cooled gas that 
is infalling toward the central galaxy. The hot wind ejected 
by a galaxy acts preferentially on the most pervasive hot 
phase, affecting in a much weaker way the cold infalling gas, 
which naturally fragments into clouds with a low covering 
factor. So, we deem that treating the infalling gas as be- 
longing to a different phase is a step forward in the physical 
description of galaxy formation, especially in those infall- 
dominated halos where a high fraction of halo gas is in the 
cold phase. As a matter of fact, the recipe by Croton et al. 
(2006) discussed above (section 15. 2p goes in the same di- 
rection, because with their assumption of A/ cosmo = Mi n ,H 
(valid for infall-dominated halos) the cooled gas is not af- 
fected by feedback. 

The cold gas is let infall to the central galaxy on a 
number nd yn of dynamical times computed at the cooling 
radius r coo i: 

fe? , (3D 



M in ,H 



M in , B = M in , H X 
Min.D = Min.H x exp 



exp 



( 2R D ) ( 1+ 2R D ) 
1 



( 2R D ){ 



Rb \ 
2Ru) 



(35) 
(36) 



Here and in the following, Rb denotes the half-mass radius 
of the bulge, while Ro is the scale radius of the disc. 

This option has a significant effect on the ability of 
quenching cooling at low redshift by AGN jets fsection l9.2|) : 
if the infalling gas has to wait for an external trigger (like 
a merger) to get into the bulge component, and from there 
to accrete onto the central BH, the feedback from the AGN 
would be activated with a significant delay with respect to 
the start of the cooling flow, while the activation is much 
quicker if the infalling gas is allowed to flow directly into 
the bulge. 

An interesting multi-phase model for cooling and infall 
has been proposed by Mailer & Bullock (2004). While part 
of the gas that resides within the cooling radius cools and 
fragments into clouds, a fair fraction of it remains at the 
same temperature and, due to the drop in density, with a 
long enough cooling time to prevent its cooling. The infall 
of the clouds to the galaxy is then followed in detail, taking 
into account the physical processes (mainly cloud collisions 
and ram-pressure drag) that make the gas loose enough ki- 
netic energy to fall into the galaxy. This process leads to 
a significant slowing down of the infall. Unfortunately, the 
Mailer & Bullock model does not take into account the effect 
of the residual pressure on the evolution of the cooling ra- 
dius. Clearly our nd yn parameter gives a poor representation 
of this complexity, and a further sophistication of the mod- 
eling of this process in MORGANA may be worth performing 
in the future. 



5.5 Galaxy super-winds and cosmological fall-back 

Whenever the gas phases of the halo component are too 
energetic to be bound to the DM halo, they are allowed to 
escape as a galaxy super-wind. 
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The hot gas is let flow away whenever its energy over- 
takes the virial one by a factor /wind: 



hot wind condition : J5hot,H > /wind-Evir 



(37) 



where E viv = (-0.5)UnM h>a /MH (see section [572)1 . The pa- 
rameter / w i n d is inserted to avoid the excessive escape of gas 
that we have noticed when it is set to unity; the results are 
rather stable when / W md ~ 2, so we adopt the best-fit value 
1.7 for the rest of the paper. Clearly this parameter should 
be fixed by a careful comparison to hydro simulations. Call- 
ing tsound = rhlc a the sound-crossing time of the halo, the 
hot wind mass and energy flows are then computed as: 



M hsw 



/wind-E'vir \ Mhot,H 



^sound 



3kT a 

j 



(38) 
(39) 



A quick computation can show that the loss of thermal 
energy by adiabatic expansion of the hot gas due to the 
mass loss should be equal to 2/3 of the energy loss: if P = 
pkT / jj, hot m p and dV = pdM then PdV = kTdA4/ fj, hot m p . 
So we set: 



-^ad ^ ^disw 



(40) 



A similar model is applied to the cold wind (with er^ = 
2_Kh/A^c,h the velocity dispersion of cold halo clouds and 
tkin = Th/o-a). The cold gas is ejected out of the halo if: 



1 2 

cold wind condition : Ku > /wind— M Ct uV^i B 
The resulting mass and energy flows are: 



M ca , 
it 



/windVd is p \ M Cj H 

&U J tkin 
/windVj; sp \ Ku 
tkin 



Ik 

^ Ji csw 



(41) 

(42) 

(43) 
(44) 



The mass ejected by the DM halo is then re-acquired 
back by it at a later time. We estimate the fall-back time as 
follows. The cold and hot gas phases escape because their 
typical velocity (kinetic or thermal) is larger than the escape 
velocity of the halo they belong to. At the end of the inte- 
gration over a time interval, we then scroll the merger tree 
forward in time and compute the time at which the parent 
DM halo has a larger escape velocity than the typical veloc- 
ity of the ejected gas. We then let this gas fall back to the 
DM halo by adding it to the cosmological infall flow (equa- 
tion [19]). However, while the large-scale structure outside a 
DM halo is clustered, galactic super-winds are emitted in 
a much more isotropic way. As a consequence, much mass 
could be ejected into voids and never fall back to the DM 
halo. We take this into account by letting only a fraction 
/back of the ejected gas fall back to the DM halo. Our results 
are remarkably insensitive of the value of this parameter; we 
use 0.5 in the following. 



6 BULGE AND DISC STRUCTURE 

For each exponential disc we record its scale radius Ro and 
its velocity Vfo at the optical radius, defined as 3.2_Rd (Per- 
sic, Salucci & Stel 1996). The half-mass radius of the disc is 
then equal to 1.6783i?D- For each bulge we record its half- 
mass radius Ru (equal to 1.35 effective radii) and its circu- 
lar velocity defined as = GMb / Rb ■ These quantities are 
sampled in the time grid defined in section [2~2l 



6.1 Discs 

The size of galaxy discs is computed following an extension 
of the model by Mo, Mao & White (1998) that takes into ac- 
count the presence of a bulge. It is assumed that the hot gas 
has the same specific angular momentum as the DM halo, 
and that this angular momentum is conserved during the 
infall. Moreover, it is assumed that the disc is exponential. 
The angular momentum of the disc is: 

poo 

Jd = 27t/ V Iot (r)X u (r)r 2 dr. (45) 
Jo 

where £d(?") = Md exp(— r / Rd) /2-kR^, is the exponential 
profile for surface density. The rotational velocity given in 
this formula contains contributions from the DM halo, bulge 
and disc: V? ot = Vff + Vb + V%. The halo contribution is 
simply Vi{r) = GMu{r)/r x (1 - fi„/fio0 

and an analo- 
gous expression is valid for the bulge, for which we assume 
a Young (1976) density profile (whose projection gives the 
observed de Vaucouleurs profile) . The disc contribution is as 
usual: 



(46) 



where y = r/2Ru and the functions contained in the equa- 
tion are the standard Bessel functions. The specific an- 
gular momentum must be equal to that of the DM halo, 
Jn/M-p = Jh/Mh- This translates into an equation for _Rd 
that must be solved iteratively, starting from the approxi- 
mate solution Ru ~ 0.71ArH (the simplest case described by 
Mo, Mao & White 1998). This computation is a bottleneck 
for the whole code, especially if disc structure is updated 
at each time-step as the profile of the hot halo gas is. To 
speed up the code, disc structure is updated each time the 
disc grows in mass by some fraction, set to 1 per cent. We 
have verified that this approximation reproduces with fair 
accuracy a disc which is growing in mass by continuous in- 
fall. Because of feedback, a disc that receives no infalling gas 
decreases in mass by ejecting a wind to the halo. This mass 
ejection presumably leads to a decrease of surface density 
with no change in the radius. We then decide to re-compute 
the radius only when the disc mass increases with respect 
to the value at the last update. 

Adiabatic contraction of the DM halo as the baryons 
settle in the centre is introduced (again following Mo, Mao 
& White 1998) by assuming that the adiabatic invariant 



B In this case the density profiles of DM and baryons are so dif- 
ferent that they need to be treated differently, so that, at variance 
with the computation of the hot halo gas profile (section 15. It and 
in agreement with Mo, Mao & White (1998), we exclude baryons 
from the computation of the DM density profile. 
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GM(r)r is constant. This implies that the DM mass within 
a radius r comes from a larger radius r; such that: 

M DM (n) + M D (r) + M B (r) = -M DM (n) (47) 

r 

(here Mdm is the unperturbed mass density profile of the 
DM halo, and only the non-baryonic DM is considered). 
Equation 1471 which must be solved together with equa- 
tions [45] and 1461 introduces a second iteration in the com- 
putation, and then a further slow-down. A solution for this 
would be to solve the equation in a 4D grid of values of the 
parameters A, c n f w , Mo/Mh and Mb/Mh, then interpolat- 
ing the solution on the table. This improvement is in project. 
In the meantime, the computation of adiabatic contraction 
is switched on only on request. We find that, at variance 
with many other authors, its introduction does not influ- 
ence strongly the results of the code, and we ascribe this 
difference to the introduction of the bulge in the compu- 
tation of the angular momentum of the disc (equation 145 [ I. 
which alone leads to more concentrated discs. 

The specific angular momentum of DM halos changes 
during their evolution, both in modulus and in direction. 
The change in the latter quantity, for instance, drives those 
precessions of discs that are commonly found in N-body sim- 
ulations. As mentioned above, PINOCCHIO can provide infor- 
mation on the angular momentum of DM halos, so it could 
be used to track its evolution. A simpler choice would be that 
of re-drawing A from its distribution (equation [9} at each 
major merger of DM halos. However, major changes in Jh 
are not handled easily within the Mo, Mao & White model 
if a disc is already formed. In fact, as the DM halo grows 
in mass the main contribution to its angular momentum is 
given by the most recently accreted mass shells, while the 
model does not take into account the internal distribution of 
angular momentum. As a consequence, the change induced 
by last accreted and not yet cooled shell would force a re- 
computation of the whole disc structure. This can lead to 
unphysical events like sudden bursts of star formation due to 
a decrease of A after a DM halo major merger. Such events 
can be avoided either by simply keeping fixed the value of 
A for each DM halo, which is our choice, or by using a more 
sophisticated algorithm for disc structure, able to consider 
the contribution to the angular momentum of the disc by 
accreting shells of gas. 

A connected sophistication lies in taking into account 
the internal distribution of angular momentum, which be- 
haves like J D (r) oc M (r) a (Warren et al. 1992; Bullock et al. 
2001), with the exponent a ranging from 1 to 1.3. This has 
important implications in the modeling of discs, as it im- 
plies that the first cooled gas has a low angular momentum, 
and then settles into a more compact, higher density disc. 
However, a straighforward implementation of this criterion 
is not easy, as it leads to a coupling of Jd to the amount 
of cooled gas and then to the disc mass. This is determined 
by the combined action of cooling and feedback, which in 
turn depends on the surface density of cold gas and then 
on Ro and Jd- This leads to nasty oscillations or numerical 
instabilities in the integration. 

This illustrates an extreme case of another issue which 
is not addressed by this code, namely that the re-heated gas 
of galaxy winds is assumed to have the same specific angular 
momentum as the halo, something that may be unrealistic in 



many cases. Clearly the distribution of angular momentum 
of gas is a topic that needs much attention. Hydro simula- 
tions are the right tool to address this issue, and yet the 
angular momentum of stellar discs is only recently showing 
some hints of numerical convergence in the biggest simula- 
tions (see, e.g., Governato et al. 2006, but see also D'Onghia 
et al. 2006). From this point of view, it is remarkable that 
the simplest assumptions on the distribution of angular mo- 
mentum give discs whose properties are not so different from 
reality. 



6.2 Bulges 

If discs contain the baryons that have retained their angu- 
lar momentum, bulges contain the baryons that have lost 
most of it. In MORGANA there are four ways to let gas lose 
angular momentum and flow to the bulge. The first is di- 
rect infall from the halo. Indeed, as suggested by Granato 
et al. (2004) and mentioned above, the first fraction of gas 
that collapses is likely to have a low angular momentum, 
and is then doomed to become a spheroid. We do not im- 
plement this idea directly, to avoid the numerical problems 
mentioned in the previous subsection. However, as explained 
in section \5A\ we allow gas to infall in an existing bulge by 
a fraction equal to the fraction of disc mass embedded in a 
bulge. 

The second mechanism is bar instability. Whenever a 
disc embedded in a DM halo has a sufficiently high sur- 
face density, it becomes unstable to bar formation. The bar 
brings a fair fraction of the disc mass into the bulge; we as- 
sume this fraction to be /bar = 0.5. For the condition for 
disc stability we use the following (Efstathiou, Lake & Ne- 
groponte 1982; Christodoulou, Shlosman & Tohline 1995; 
Mo, Mao & White 1998): 

The value for eiimit ranges from 1.1 for stellar discs to 0.9 
for gaseous disc; as this effect is important at high red- 
shift, when discs are mostly gaseous, we use the lower value. 
This condition is checked at each integration time-step, and 
whenever it is not satisfied the integration is interrupted. At 
this point a fraction /b ar of the disc mass (that we assume to 
be 0.5) is given to the bulge, and the integration is started 
again; in other words, bar instability is implemented as an 
external flow. No starburst is explicitly connected with this 
event, but the presence of gas in the bulge gives rise to a 
stronger burst of star formation (see section r7.3[) . The re- 
sults are not strongly sensitive to the parameter /bar, as the 
unstable discs typically acquire so much mass that they un- 
dergo a series of consecutive bar instabilities; changing the 
value of /bar influences a bit the pattern of repeated star- 
bursts but does not influence much the final mass in the 
bulge. 

The third mechanism is of course the merging of galax- 
ies. The implementation of mergers has already been de- 
scribed in section 14.31 It is worth repeating here that the 
uncertainty in the fate of gas in minor mergers (we put the 
whole satellite in the bulge, other authors make different 
choices) does not influence much the results. 

In both cases of disc instability and mergers the radius 
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Rb of the bulge formed by the merger of objects 1 and 2 is 
computed as suggested by Cole et al. (2000): 

(Mi + M 2 ) 2 _ Ml | Ml | / MxM 2 
Tie ^i c i?i + i?2 

Here i?i and i?2 are the half-mass radii of the merging galax- 
ies, computed assuming (as done above) an exponential pro- 
file for the disc and a Young (1976) profile for the bulge. For 
the parameter //ewe use the values of 2 for mergers and 4 
for disc instabilities, as suggested by Cole et al. (2000). 

The fourth mechanism to form a bulge is through feed- 
back. A necessary condition for a gaseous disc not to transfer 
angular momentum outwards is that the dissipative cold gas 
has a low enough velocity dispersion, a condition which is 
satisfied in observed nearby discs. However, the kinetic en- 
ergy of the cold phase is determined by feedback from SNe. 
We anticipate (see section 17. ip that for discs with a high 
enough surface density of cold gas the velocity dispersion of 
cold clouds can be much higher. This process is observed to 
be at play in high redshift galaxies (Genzel et al. 2006) and 
implies a loss of angular momentum and a corresponding 
thickening of the disc into a bulgy object. We model this 
event by simply forcing a bar instability (i.e. transferring 
a fraction /bar of mass to the bulge) whenever the surface 
density of cold gas E co id,D overtakes a limiting value: 

Scold, D > Slimit (50) 

This mechanism is used in Fontanot et al. (2006a) to enhance 
the number of bright quasars, but is not used in the results 
presented in this paper. 




Log H c „ (po) 

Figure 5. Feedback regimes as a function of surface mass den- 
sity and vertical scale-length. Green stars denote systems in the 
adiabatic blow-out regimes (the encircled points highlight the un- 
stable systems; see M04), red dots denote systems in the adiabatic 
confinement regime, blue crosses denote systems in the PDS con- 
finement regime. The blue continuous lines separate the regions 
with different feedback regimes. The green and red shaded areas 
give the typical regions occupied by galaxy discs and bulges. This 
is an adapted version of figure 3a of M04. 



7 STAR FORMATION AND FEEDBACK 

The treatment of feedback in the present version of MOR- 
GANA follows the multi-phase model of M04. As anticipated 
in section 12.31 instead of implementing a full multi-phase 
model of the ISM (which would be straightforward as far 
as the coding is concerned, but would lead to a number of 
numerical problems) we prefer to insert the main results of 
the M04 model as simple recipes, so as to simplify the nu- 
merical integration and to gain a more immediate grasp on 
the physical processes inserted. Anyway, we stress that the 
insertion of a physically motivated model for feedback in 
place of a set of phenomenological recipes (as done in most 
other galaxy formation models) is one of the most important 
features of MORGANA. 

M04 studied the evolution of the ISM under the follow- 
ing assumptions: 

• the ISM is composed by two phases, a hot and a cold 
one, in thermal pressure equilibrium; 

• collapsing and star-forming clouds arise from the cold 
clouds by kinetic aggregations; 

• type II SNe exploding within a star-forming cloud give 
rise to a single super-bubble; 

• the super-bubble propagates in the most pervasive hot 
phase; 

• the super-bubble expands until either it is stopped by 
the external pressure or it blows out of the structure. 

Four possible self-regulated feedback regimes follow natu- 
rally from this setting, depending on whether the super- 



bubbles stop by pressure confinement or blow-out, before or 
after the internal hot gas has started cooling (so as to form 
a pressure-driven snowploughs, hereafter PDS). 

In M04 it was shown that the dynamics of feedback de- 
pends mainly on the total surface density of the disc Ed and 
vertical scale-height H e ff of the system. In realistic systems 
feedback can take place in the regimes where super-bubbles 
blow-out or are confined in the adiabatic stage. Figure [5] 
shows that the adiabatic blow-out and confinement regimes 
take place in two different regions of the Y,o-H e ff space, sep- 
arated by the relation: 

ED - 8 (»fe)" 08 M0PC_2 (51) 

The numerical constants in this relation depend on the (un- 
certain) values of the parameters used in the model, so they 
are to be considered as indicative. 

Galaxy discs and bulges tend to stay respectively in the 
adiabatic blow-out and confinement regimes, so at a basic 
level there is no much doubt on how to apply the feedback 
regimes within a galaxy formation code: discs are in the 
adiabatic blow-out regime, bulges are in the adiabatic con- 
finement regime. However, the central region of a disc can be 
in the adiabatic confinement regime for at least two reasons: 
(i) it is embedded in a bulge, so that its ISM is pressurized 
by the bulge hot phase; (ii) its surface density is high enough 
to cross the limit of equation [51] for adiabatic confinement. 
This happens at a typical surface density of Ed ~ 100 — 300 
M© pc~ 2 . Point (i) has been used to justify the direct infall 
of gas from the halo to the bulge (section [5]4}. Point (ii) has 
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been mentioned in section [6721 to argue for feedback-induced 
loss of angular momentum, and the effects of such feedback 
have been modeled by forcing bar formation when the sur- 
face density of cold gas is high (equation I50[) . We show more 
details and results on this mechanisms in Fontanot et al. 
(2006a). 

To fully implement the M04 feedback recipes in MOR- 
GANA it is necessary to assess the role of the velocity dis- 
persion of cold clouds, a quantity which is left as a free 
parameter in that paper. This is done in the following sub- 
section. After that we explain how we implement feedback 
in the case of thin or thick systems. 



7.1 Kinetic feedback and the velocity dispersion 
of cold clouds 

We do not include here a full self-consistent treatment of the 
kinetic energy of cold clouds, but try to take into account its 
effect as follows. The velocity dispersion cr co id = y/2K/McxM 
of cold clouds in the ISM of a star-forming galaxy is deter- 
mined as a first approximation (i.e. neglecting all mass flows) 
by the equilibrium between the injection of kinetic energy 
by SNe and the dissipation by turbulence: 



K = Ksn — K dB 

The injection of kinetic energy is: 
• Met p 

ASN = e k-rj -C/SN 

-M*.SN 



(52) 



(53) 



where Esn is the energy of a SN, M* t sN is the mass of newly 
formed stars per SN and is its fraction given to the cold 
phase as kinetic energy. The rate of energy dissipation is 
observed in hydro simulation to be of the order (Mac Low 
2003): 



Ki 



1 M co ld <T cold 



(54) 



2 ^drive 

where Ldrive is the scale at which turbulence is driven, sug- 
gested to be twice the diameter of the typical size of the 
super-bubble. To determine o- C oid we assume that an equi- 
librium configuration is quickly reached, so that K = 0. We 
then obtain: 

( u V 1/3 

o-cold = cr — - — (55) 

V 1 Gyr / 

where <ro = [2Z/d r ivc£fc-EsN/M*,sN(l Gyr)] 1 / 3 and the time- 
scale for star formation is t* = M co \&/ M B f. In a thin system 
like the Milky Way the fraction of (thermal and kinetic) 
energy injected into the ISM is predicted by M04 to be about 
5 per cent, so et, will be of order of 0.01, while the driving 
scale should be of order of the vertical scale-length of the 
system, ~ 100 pc. The numerical value of ctq is then: 



g o = 9.3 



100 pc 

1/3 



Vo.oi/ 



1/3 



1/3 



(56) 



E ' 



KL 



SN 



120 M 



-1/3 



km s 



The predicted value of cr co id for a Milky Way-like galaxy 
(with t* ~ 2 Gyr) is then ~ 7 km s~ , in remarkable agree- 
ment with the observed value which ranges between 6 and 9 



km s _1 . Much higher values of the velocity dispersion are ex- 
pected in thick systems, where et is much higher (all the SN 
energy is retained by the system) and t+ much lower. This is 
in agreement with the sparse data available (see, e.g., Dib, 
Bell & Burkert 2006 for local galaxies, and Genzel et al. 2006 
for galaxies at z ~ 2). As a consequence, we neglect kinetic 
feedback in discs but take it into account in bulges. Given 
the great uncertainty in the driving scale of turbulence in 
mergers, the normalization constant <to will be left as a free 
parameter. 

We stress that equation [56] is valid under the assump- 
tion that turbulence is driven by SNe and not by gravita- 
tionally induced motions like differential rotation or tidal 
disturbances in mergers. In the case of our Galaxy, Mac 
Low (2003) argues from energetic arguments that SNe are 
the most likely drivers of turbulence, and the good agree- 
ment of the predicted value of <r C oid given above with the 
observed one suggests that they are at least a significant con- 
tributor. In thick systems, like gas-rich spheroids or mergers, 
gravitationally-induccd turbulence will be very important. 
However, the most relevant consequence of a high value of 
0"coid for our purposes is the massive removal of cold gas 
during episodes of strong star formation. It is clear that if a 
condition cr co id > Vb is reached in a gas-rich bulge, then the 
driver of turbulence must be star formation and not gravity, 
while in the case of a low value of o" co id no mass is removed. 
So, in this case the assumption will work in the range where 
its results are most influential. 



7.2 Discs as thin systems 

In thin systems, super-bubbles blow out of the disc soon 
after they form and while they are still in the adiabatic stage; 
their porosity remains low. The fraction of the SN energy 
that is injected into the ISM is limited to a few per cent of 
the total budget, while most energy is directly injected into 
the halo through a tenuous, metal-rich hot wind. The ISM 
self-regulates to a configuration that is very similar to that 
of the Milky Way (M04). The star formation regulates to 
a level similar to that found in nearby galaxies by Schmidt 
(1959) and quantified more recently by Kennicutt (1998). 

In M04 the resulting star-formation timescales were in- 
correctly scaled with the infall time and the total surface 
density Ed of the system. A more careful analysis reveals 
that, for a fixed set of parameters and at fixed H e g, de- 
pends mostly on the cold gas surface density E co id,D- This is 
expected, as the star-formation timescale is determined by 
the intrinsic properties of the ISM, not by the rate at which 
mass is acquired by it or by the amount of stars present. 
A better fit of the results of M04 in the adiabatic blow-out 
regime (for the standard set of parameters defined in that 
paper) gives: 



U = 1 



1 M 



pc" 



100 pc 



Gyr 



(57) 



In the simple case of a constant value of H e fs for all 
galaxies, the Schmidt-Kennicutt law is recovered with a very 
similar normalization (however, observations refer to inte- 
grated quantities, this prediction refers to a section of the 
disc). But the assumption of a constant H c g is rather ar- 
tificial, as the vertical scale-length of the cold gas is de- 
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termined by the vertical gravitational equilibrium as H e fi = 
o-^ oldiD /7rGED, so a constant value would imply a tuned vari- 
ation of Ed and <x C oid,D- The other simple possibility is to 
assume a constant value for cr co id,D' This leads to the predic- 
tion t* oc E ra } d ° D /co?d> where / co id,D is the fraction of cold 
gas in the disc. This relation is apparently steeper than the 
observed Schmidt law; however, E co id,D and / co id,D are usu- 
ally correlated in galaxies, in that lower E co id,D discs, having 
lower star-formation rates, retain a greater fraction of their 
gas. For instance, with our reference choice of parameters 
(section I12[) we obtain / co id,D oc E° D * d D) which implies the 
much shallower relation t* oc D . 

The dependence of the velocity dispersion on the star- 
formation timescale, equation 1551 is then inserted, through 
the definition of H c g, into equation 1571 and the resulting 
relation is normalized so as to lie on the average Schmidt- 
Kennicutt relation at E co id.D = 13 Mq pc~ 2 (the gas density 
of the Milky Way, obtained by assuming Ro — 3.5 kpc and 
M co i dj D = 5 x 10 9 M Q , see Cox 2000). The resulting relation 
for the star-formation timescale in discs is: 

**- D = 9 - 1 l lMopc-^ l-bl-J Gyr (58) 

This relation is again apparently steeper than the Schmidt- 
Kennicutt one, but when the correlation between E co id,D 
and /cow is taken into account we obtain t+,D oc E^j 3 ^ or 
E s fr,D oc Ej ^d,D) m excellent agreement with the observed 
Schmidt law. 

Other complications can be introduced in this picture: if 
the driving length Ld r i vc is equal to twice the diameter of the 
blowing-out super-bubbles, it will scale with H a g. Moreover, 
the efficiency of feedback is known in M04 to increase slightly 
with Ed (and then with E co id,D)- Introducing these depen- 
dences we obtain a scaling of t± with E co id,D and / co id,D 
similar to equation [55] with slightly different exponents. We 
have verified that these somewhat different descriptions of 
feedback do not give significantly different results, so we keep 
equation [55] as a reference. 

According to M04, in the adiabatic blow-out regime the 
mass ejection rate due to the action of super-bubbles is mod- 
est, while much gas belonging to the hot phase leaks out of 
the halo simply because it is too hot to be confined by the 
disc. This unavoidable mass flow, observed as a hot corona 
surrounding the star-forming discs (see, e.g., Fraternali et 
al. 2002), amounts to a mass ejection rate very similar to 
the star-formation rate. We then set the disc hot wind term 
A^hw.D as equal to the star formation rate M s f,D- This is at 
variance with most other galaxy formation models, where 
the rate of gas re- heating, equivalent to our Mh w ,D , is related 
to the star-formation rate through a /3 parameter, assumed 
to scale as a power of the disc velocity. Our assumption is 
equivalent to /3 = 1; the higher ejection efficiency in small 
discs will then be due to galaxy super-winds triggered by the 
energy injection from the galaxy. Finally, we neglect any cold 
wind term, due to the results of M04 and to the predicted 
low level of turbulence in discs (section l7.ip . 

Using the instantaneous recycling approximation (see 
also section [5}, we let a fraction / rcst of the mass of newly 
formed stars to be restored to the ISM. Summing up, the 
star-formation, restoration and wind mass flows in discs are: 



M a f } T) = M co id,D/i*,D 

Afrs.D = /rest A{f s f,D 

A^hw.D = A-^sf.D 

M CW , D = 

The hot wind mass flow does not carry much thermal 
energy: assuming a typical temperature of 10 6 K (M04) and 
for M*,sn = 120 M and 10 51 erg per SN, the fraction of 
SN energy carried away by this wind is ~ 0.05. However, 
most SN energy is blown out to the halo by the blowing- 
out super-bubbles, and this energy mixes with the thermal 
energy of the hot outflowing wind. We then assume that a 
fraction /th,D of a 10 51 erg SN is carried by the hot wind gas; 
this fraction is estimated to be ~ 0.8 by Monaco (2004b). 
However, this number is highly uncertain because (i) the loss 
of SN energy within the star-forming cloud could be higher, 
as many authors working on the physics of ISM have often 
claimed, (ii) the SN energy could be well in excess of 10 51 
erg, (iii) some energy could be lost in the interaction with 
the hot halo gas. Moreover, there is a degeneracy between 
this parameter and M*,sn, which depends on the IMF. It is 
then wise to leave /th,D as a free parameter; it is however 
remarkable that its best-fit value turns out to be 0.5, not far 
from the value suggested by Monaco (2004b). 

To the energy of type II SNe we add a contribution from 
type la SNe as one per year per 10 12 Mq of stars. This is 
done in order to test the energetic effect of these SNe. In 
the case of thin systems with a continuous star formation, 
and given the uncertainty on /th,D, the contribution from 
type la SNe is almost irrelevant. As suggested for instance 
by Recchi et al. (2002), type la SNe may have a higher effi- 
ciency of energy injection, as they don't explode in the dense 
molecular clouds. As a consequence, one might want to have 
a different efficiency /th,ia for the two SN types. However, 
given the modest importance of these objects we prefer not 
to add a further parameter. 

The contribution of the disc to the hot wind energy 
term is then: 

E hw ,nL = VdSsn + 1 SN yr-i-^f ) (60) 

ldlsc \M*,sn 10 12 M J 

The outflowing hot gas will interact with the halo cold phase 
in a way that is difficult to model. In order to be able to 
constrain the effect of this process, we allow for a fraction 
of the thermal energy to be given to the halo cold phase: 

-Kw.H = / kin i?hw,H ,. (61) 

I disc I disc 

For simplicity we do not subtract this energy from the hot 
wind budget. We stress that this term is not due to a con- 
tribution of cold wind but to the interaction of the hot wind 
with the cold halo phase. 

Star formation in the external parts of galaxy discs is 
known to be truncated (see, e.g., Kennicutt 1989), and this 
is likely due to the differential rotation, even though also 
the thermal transition from warm to cold gas may play a 
crucial role, as suggested by Shaye (2004). This author also 
suggests that the threshold for star formation may be sim- 
ply expressed in terms of a threshold for gas surface density, 
which scales with the disc and ISM properties and, very im- 
portantly, with the ionizing UV background. We implement 
this threshold E t hr for star formation simply by increasing 
the star-formation timescale t*,D by the inverse of the frac- 
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tion of disc mass for which E co id,D > Ethr- Clearly, a better 
but more complicated choice would be to divide the disc into 
three zones, corresponding to the bulge, the star-forming 
disc and the gas beyond the star-formation threshold; this 
is left to future work. We have also tried to modulate this 
threshold with redshift, following the increase by an order 
of magnitude of the ionizing background from z = to 1 
(Bianchi, Cristiani & Kim 2001). However, the dependence 
of the threshold on the ionizing background, is so weak that 
this modulation does not influence the results appreciably. 



7.3 Bulges as thick systems 

In the model by M04 the gravitational perturbations to the 
ISM are neglected. This assumption is questionable for spiral 
discs, but the coincidental similarity between the timescale 
for kinetic aggregations and the disc timescale, that drives 
the sweeping by spiral arms, makes the kinetic aggregation 
mechanism a good substitute for spiral arms in creating mas- 
sive clouds, at least at an order-of-magnitude level. Things 
change in the case of thick systems like bulges, where most 
energy is efficiently pumped into the hot phase so that the 
ISM is much more pressurized, the star-forming clouds are 
much smaller and denser, and kinetic aggregations are dis- 
favored. This leads to the prediction of a lower level of star 
formation compared with the Schmidt-Kennicutt law (figure 
7 of M04). However, tidal disturbances will be much stronger 
in a bulgy object, at least when it is formed through a disc 
instability or a merger, so that the formation of collapsing 
clouds will not be determined by kinetic aggregations. More- 
over, in a major merger the pressurization due to the onset 
of the adiabatic confinement regime will cause a quick drop 
of the Jeans mass, and this will make most clouds present 
in the merging discs collapse and form stars in a short time. 
These transient effects have not been properly modeled by 
M04, so we prefer to determine the star-formation timescale 
directly from the Schmidt law: 
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(62) 



The resulting star-formation and restoration rates in bulges 
are: 



Afrf.B 
M Ia B 



M c ,bA*,b 

/rest M B { z B 



(63) 



As for the thin system case, the rate at which hot gas 
flows to the halo is similar to the star formation rate, so 
that we retain this prediction. While in thin systems most 
energy (though not most mass) is injected to the halo mainly 
by blowing-out super-bubbles, in thick systems the energy 
is ejected mainly through this hot wind. In the M04 model 
the typical temperature of the hot phase is found to be of 
~ 10 r K, higher by an order of magnitude than the thin 
system case; the fraction of SN energy carried away by the 
hot wind is then /th,B — 0.5, slightly lower than the thin 
system case. Moreover, while a disc is unable to confine a 
gas phase with a temperature as high as ~ 10 6 K, a massive 
bulge can gravitationally confine its hot phase, whose tem- 
perature corresponds roughly to the virial temperature of a 
bulge with Vb ~ Vh ot = 300 km s _1 . We then limit the hot 
wind as follows: 



2 

hot 



(64) 



where of course (3 = if Vb > Vhot- The resulting hot wind 
rate is: 



Mhw.B L = /3M sf , B 

I therm 

while the hot wind energy is: 
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bulge 



(65) 
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As for the thin systems, we allow some energy to accelerate 
the cold halo phase: 



K, 



w,H 



bulge, hw 



/kin -Ehw.H 



bulge 



(67) 



Because in thick systems all the energy of SNe is in- 
jected to the ISM and star-formation timescales are much 
shorter (in virtue of the increased surface density), the in- 
jection of kinetic energy leads to velocity dispersions of the 
cold phase much larger than the ~7 km s _1 found in discs. 
In this case some cold clouds may get enough kinetic energy 
to be able to leave the potential well of the bulge. The prob- 
ability that a cold cloud is unbound, i.e. it has a velocity 
larger than the escape velocity of the bulge, \P2Vb-, is com- 
puted under the hypothesis of a Maxwellian distribution of 
velocities with rms Ucoid b: 



funb = erfc 
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(68) 



The average velocity of the unbound clouds is 

2 



Wunb — 0"cold,B 



1+ , 

4 \ \ CT co i d ,B 



exp 



2?r 



(69) 



The velocity of the clouds after being ejected out of the bulge 
will be: 



2V£ 



(70) 



This process will generate a wind flux roughly equal to 
the unbound mass, P un bAfc,B, divided by the crossing time 
flWtWb- This outflow is naturally identified with a cold 
wind, adding mass to the cold halo phase. The cold wind 
mass flow is: 



I bulge, cw 



M c , B Pur 
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(71) 
(72) 



However, the fate of the cold clouds ejected by a bulge may 
be different. The exit from the pressurized bulge to the halo 
would lead to an expansion of those clouds, making them 
much more sensitive to the bow shock generated in their 
interaction with the hot halo gas and to thermal evapora- 
tion. As a result, the cold clouds ejected by the bulge could 
be heated and mix with the hot phase. This is analogous 
(though not in the details) to what happens when kinetic 
feedback is implemented in SPH simulations: a cold parti- 
cle neighbouring a star-forming region is given some kinetic 
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energy, but this energy is quickly thermalized by the inter- 
action with the other particles. If this is the case, the cold 
and hot wind flows should be given as follows: 



M h 



I kin 

M hw , B 

Mcw.B 



Mc.bPui 



«unb 



MwB L + Afhw.B , . 

I therm I kin 



= 



(73) 
(74) 
(75) 



From the energetic point of view, the kinetic energy of this 
outflow will be typically not enough to heat the gas to a high 
temperature, so the heating will be done with the same en- 
ergy budget of equation [67] in other words, no energy is 
added to the -E n w,H term. We have implemented both pos- 
sibilities, and the choice between the two has been left free. 
The results presented here are obtained giving the outflow- 
ing gas to the hot phase. 



8 METALS 

The evolution of metals is given by the equations reported 
in table 1. In this set of equations most metal flows are 
obtained from their related mass flow as follows: 



Mff ow 



AL 



-M flo 



(76) 



where M sourcc and M source are the gas and metal masses of 
the source phase. The cosmological infalling gas is assumed 



to have a metallicity Z pT , 



10" b due to pre-enrichment 



by sources that are below the mass resolution limit (from 
popIII stars to very small primeval galaxies). These flows 
take into account the transfer of metals among phases and 
components, but not their production. 

We assume that (i) metals are produced by newly 
formed stars in the instantaneous recycling approximation, 
i.e. their production follows instantaneously star formation; 
(ii) the new metals are instantaneously mixed with the ISM. 
Newly metals are spread into the ISM mainly by SNe, so 
analogously to the energy they are likely to be selectively 
ejected to the halo. This is clearly true in the blowing- 
out thin systems, but even in thick systems metals are first 
mixed with the hot gas that escapes the halo at a rate equal 
to the star-formation rate. To model this effect without us- 
ing explicitly a multi-phase description of the ISM, a fraction 
fzej of the new metals is ejected directly to the halo through 
hot winds. In the case of bulges, this fraction is multiplied 
by the /3 factor of equation l64l to take into account the abil- 
ity of massive bulges to retain the outflowing hot gas. We 
have: 



Mi 



jy. y (Mf,i 



M sf , D 



(77) 



where Y is the fraction of mass in newly produced metals 
per generation of stars. This term is added to the M hW]H 
metal flow; satellite galaxies will inject their metals to the 
halo component of the main DM halo they belong to through 
the satellite metal flows. The other metals will be given to 
the ISM of the component they belong to: 



A value of fzej = 0.5 will be used in the following. 

9 AGN ACTIVITY 

9.1 Accretion onto black holes 

The modeling of black hole accretion in MORGANA has al- 
ready been described, though in a slightly different way, by 
Monaco & Fontanot (2005), and will be briefly summarized 
here. 

A seed black hole of mass M Bee & is put at the centre of 
each DM halo. These black holes may be generated by the 
collapse of the first stars (e.g., Volonteri, Haardt & Madau 
2003). Seed masses should be of order of tens to hundreds 
Mq ; however, they start growing very soon during the early 
evolution of baryons in DM halos. This happens at times 
and for DM halo masses that are not sampled in the typical 
runs used in galaxy formation. We then use a higher value 
for the seed mass, M aee d = 1000 Mq; the results are quite 
stable for reasonable variations of this parameter. 

The accretion of gas onto the black holes is possible only 
if this gas has lost nearly all of its angular momentum. The 
first step in this loss is the same that leads to the forma- 
tion of bulges; we then base our computation of black hole 
accretion on the cold bulge gas. As the amount of accreted 
gas is small, we do not remove the accreted mass from the 
matter budget. In other words, the mass in black holes does 
not obey a mass conservation constraint as the mass of all 
the other components. 

The gas is assumed to lose angular momentum at a rate 
proportional to the star-formation rate. This is justified by 
the radiation drag mechanism proposed by Umemura (2001) 
and used by Granato et al. (2004) . However, a connection be- 
tween loss of angular momentum and star formation likely 
has a more general validity, as many mechanisms able to 
cause a loss of angular momentum (turbulence, kinetic ag- 
gregations etc.) are directly or indirectly driven by massive 
stars and SNe. The simplest way of modeling the loss of 
angular momentum is then: 



= /lowjM SIj B 



(79) 



More generally, the accreted mass could scale as a power law 
of the star-formation rate; for instance, angular momentum 
loss driven by encounters would likely scale with the square 
of the driving force. This is described in detail in Fontanot 
et al. (2006a); we describe here only the simplest choice. 

We follow Granato et al. (2004) by assuming that this 
gas does not flow directly onto the black hole but settles 
on a reservoir of low-angular momentum gas, whose mass is 
Mresv From this reservoir the gas accretes onto the black 
hole at a rate regulated by the viscous time-scale of the 
accretion disc, modeled by Granato et al. (2004) as: 
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(80) 



,M BH , 

where the constant fc accr is suggested by the authors to take 
(Tb — 0.65Vb is the ID velocity dispersion 



a valued ~0.001, 
of the bulge and Mbh is of course the BH mass. 



M y 1, B 

Myi D 



(l-/3/zcj)yM sf ,B 

(1 - /zcjjVA^f.D , 



7 Because of a misprint, the value is indicated to be 10 4 in their 
(78) paper, while the correct value is an order of magnitude larger. 
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Accretion is always limited by the Eddington-Salpeter 
rate, whose time-scale (for a radiation efficiency of 0.1) is 
ted = 0.04 Gyr. The evolution of the black hole-reservoir 
system is then: 



Mbh = min 

Aircsv = Afl ow j 



Mbh 



(81) 



In Monaco & Fontanot (2005), the accretion rate onto the 
black hole was modeled simply as M vesv /t^d- We prefer equa- 
tion [80] in virtue of its sounder physical motivation. 

9.2 Feedback from jets 

Jets coming from radio-loud AGNs are though to be one of 
the most promising mechanisms to stop the cooling flows 
in galaxy clusters (McNamara et al. 2005; Voit & Don- 
ahue 2005). Besides, a successful reproduction of the high- 
luminosity cutoff of the luminosity function of galaxies re- 
quires a proper modeling of this kind of feedback (Ben- 
son et al. 2003; Bower et al. 2006; Croton et al. 2006). 
The efficiency of radiation of AGNs is known to decrease 
when the accretion rate in units of the Eddington rate 
m = A/BHted/MBH is lower than ~ 10 -2 ; these slowly ac- 
creting black holes however radiate very efficiently in the 
radio (see the discussion in Merloni, Heinz & Di Matteo 
2003). It is then reasonable to assume that the efficiency of 
energy emission in jets is 0.1 if m is small, and lower for 
higher accretion rates. As ~10 per cent of bright QSOs are 
radio-loud, we estimate the efficiency of jet emission as 0.01 
in the case of rh > 0.01. 

As suggested by Croton et al. (2006), the efficiency with 
which this energy heats the hot halo gas component should 
scale with the hot gas temperature to the power 3/2. In 
order to keep the expression simple and avoid possible un- 
wanted numerical instabilities in the integration, we scale 
the efficiency with the circular velocity of the halo: 



/jet — /jot,o y-_ 



.1000 km s" 1 . 
The energy injected into the halo hot gas is then: 



/jctO.lMBHC 2 



j cts \ /jetO.OlMj 



BHC 
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m < 0.01 
m > 0.01 



(82) 



(83) 



This mechanism benefits much by the direct infall of gas to 
the bulge, described in section [6721 but, as demonstrated in 
appendix [Bj works only for particle masses not higher than 
10 9 Mq. 



9.3 Quasar-triggered galaxy winds 

A bright quasar shining into a star-forming bulge can inject 
a great amount of energy into the ISM, leading to a massive 
removal of cold gas. This mechanism has been described in 
detail by Monaco & Fontanot (2005). These winds have a 
modest effect on the formation of a galaxy, but influence the 
formation and evolution of AGNs. This is described in full 
detail in Fontanot et al. (2006a), where models with quasar- 
triggered winds are presented. We do not use these winds 
in the results presented in this paper, as the introduction 
of winds influences deeply the number of bright quasars but 
only modestly their host galaxies. 



10 PARAMETERS 

Table[4]gives a complete list of the parameters of the model, 
with the value used to compute the results given in this 
paper. To these the cosmological parameters (Qo, Oa, Sib , 
Ho and as', see the Introduction) should be added; these 
are now fixed with a good accuracy, with the exception of 
as whose value influences strongly the number of galactic 
DM halos at high redshift. The number of parameters is 
high, and this reflects both the number of physical processes 
that are included in the galaxy formation code and the level 
of uncertainty in many of these processes. Anyway, galaxy 
formation is a problem of complexity, and there is no way to 
reduce the number of parameters other than hiding them by 
fixing them to some value. On the other hand, the number of 
observables that can be used to constrain these parameters 
is very high, so their values can be fixed, with the exception 
of a few degeneracies. The most obvious one in this context 
is the degeneracy between the energy of a SN, the star mass 
per SN M ata r,sN and the thermal feedback efficiencies /th,B 
and /th,D feauations l60l and l67[) . For this reason we only vary 
the efficiencies, leaving M*,sn fixed and not even including 
the SN energy as a parameter. 

In practice, many of these parameters are fixed inde- 
pendently by the results of N-body simulations (like the pa- 
rameters of merging) or by the knowledge of the IMF, and 
their variation within the known uncertainties does not influ- 
ence the results strongly. The remaining parameters, labeled 
with a mark in the table, are then of primary importance. 
Their value is fixed by comparing with a set of basic data 
(section |12[) . we give here a very brief list of the observable 
which is most useful to fix each parameter. 

/scatter: it regulates the amount of halo stars, and is 
fixed by requiring a fraction of at least ~ 20 per cent of halo 
stars in groups and clusters (Mh > 10 14 Mq). 

Jiqucnch and Udyn- they regulate the star formation den- 
sity, especially at high redshift; their value is fixed by repro- 
ducing the very uncertain star formation density at z > 4. 
They are nearly degenerate and their value is especially sen- 
sitive to mass resolution. 

7 P : it regulates star formation, especially at low redshift, 
but its influence is modest when it is varied within the range 
suggested by observations. 

/wind: by regulating the energy level at which a galactic 
super-wind is triggered, it influences the relation between 
stellar and DM mass. 

heat cold gas: its switching on gives an effect similar 
to an increase of n qucllcn or nd yn by 1. 

close hole: when it is on, it increases cooling by a sig- 
nificant amount which depends on y p . All the parameters of 
the hot halo gas and stellar feedback should be re-tuned if 
this option is active. 

infall on bulge: it influences mostly the jet feedback, 
so it should be switched on if an efficient quenching of the 
cooling flow is wanted. 

Y: it regulates the amount of metals produced by stars, 
and is fixed by requiring a good match of galaxy metallici- 
ties. 

/zcj: it regulates the fraction of metals injected in the 
halo component, and is fixed by requiring a good match of 
the metal enrichment of the hot halo gas. 

/th.D: it regulates the stellar mass of small galaxies, and 



© RAS, MNRAS 000,[TH35] 



MORGANA 23 



Name 


reference 


Comment 


Constraint 


Equation/ 




model 






Section 






lilt." I gel S 






f, 

J n mm 


0.2 


Tnnir^T - vyipvctpt" rnnnitifin frn* T)A/f nalrtQ 

111<3j 1 VJl lllL_ltiL_l \_ U 1 1 U_ 1 U 1 U 1 1 1U1 -L/lvl HCVlUo 


N-hnd v 


en [TOl 

Cq. | J.OI 


J gmm 


3 


i n m hit nwirov^r mi ic n 1 ii it i |j"ir" i v ' i 1 • i l ."^ 
lIlcLJVJl 1I1L. I liL_ I V_.U11U1 L1U11 1U1 g CLI el A. 1 L o 


NT Kr_H"v 
in uuu y 


on rm 

L M- 1 1 _J 


fir 


2.0 /4.0 


hi 1 1 cfp TfivTyiatinn in itiot'O'ot'q /nicr mc.tnhili'hpc: 

kJ LllfciL. 1U1111<_LL1U11 111 llll_.ltiL_li_>/ UloV_. lllo UOiUlll 11L5 


NT-hod v 
i n ij uu y 


eq. 1491 


f 

J scatter 


n i 


fraction of stars scattered at a galaxy major merger 


in - uuuy 


«ortlA 11 

_eeLi^_. oi 






hello component 








1.15 


r"\ f 1 1 i,'t WM". 1 C ITinPY t \ \ trip Vi ("it era t; 
LJU1 y LI W IJ 1\_ lllLI^.A Wl Lilt. 11U L tiCl.__ 


.". 1 1QO 

U l_y__ L_ 1 V . 


cq. 1 J.OI 


f i i 

7 shock 


1.2 


Qnrir'L" ViPC-tincr fnftnr 

o lll_JV_.1V llV_.<__Llllti 1<XL_. LUl 


1 1 uu v_i y 


ens rsrn rm 

l_.V_[__ . | -_■ \J \ |_- -L | 


life! a. L Hdo 


YES 


ci.fifpn Ti^.T nn^itirifr -^/".l _"i n :-. 1 itqg c. t ty. o vy"i arcrarG 
awiiitii iui iiL_<ALinti l. vjili iieiivJ gets a l iiicljui iiit.1 a 


TV hnrh,' 

1N UUU V 


i&r\ 1^ 91 

__t.v_.L. |U ■ Zj | 


^■quench 


1 n 

± .u 


no. of crossing times for cjuenching cooling 


free 




t~~\ Ota Vi n 1 o 


YES 


suntrn trw plriQincr t hr> rv^ril i n cr n_~_lp 

□ Wl LrV_.ll 1U1 V_.lVJ__ lllti Lilt. L-VJUllllti 11VJ1L. 


free 


en f2Ql 


inf 3.1 1 on bulg© 


YES 


switch for allowing infall on the bulge 


tree 


pn« rr^i ft^i 

L_C_|b. |00| IQtJl 


^Myn 


i n 


no. of dynamical times for infall 


free 


pn e prn -TTTi 


f . , 
J wind 


1. / 


energy factor to trigger a super-wind 


free 


„^, |Q7| 1/111 

cqs. \o 1 1 rill 


fx. i 
J back 


5 


1 7' o r*i~ ~\ r"i vi vi T oiinoi" "vnriTi/H m aoo fhof t_. 1 1 o r~\ /""I/" 

11 O/L. L1U11 Ul alipLI WlilU Ill0.____ LIlO/L Iclllc. fJclv_._V 




__L_v_.L. |U ■ Ul 






disc structure 






e limit 


0.9 


limit tVvf har iriotn hilit"vr 

111111 L 1U1 VJ Cl'l lllOLCLUlllLy 


"NT-hod v 

i n u u v_i y 


en P^l 


fx. 

J bar 


0.5 


TVPrT'tir.n riT fwar tnnt otipq tr. h 1 1 1 cr> 

11 CLL. LlVJll VJl V_11__V_. LlldjL tULj L W ULlltit. 




qprt Ifi 91 

__L_V_.L. | U ■ __j | 


CLVJL laUallL i_> LJ11 IsX • 


NO 


msritr'n frw ania nnf ir /-ri vi 1 7" s r't i r~i vi 

o W ILl-ll 11.' 1 CO llil 1 Kl I H. V_. VJ11L1 cXV_. LlVJll 


free 
ree 


en [471 

cq. 1 1 






st&rs 3.nd met 3.1s 






iW *,SN 


120 

_L __ \ ) 1V1 f i j 


Q+ollm - maQG iinr ^ l\l 

__ LC>llcXl llldoo (_/ JIN 


IMF 


ens [fifTl [?T7l 

V_.L.£o, | \.}\ J I |L> 1 I 


J rest 


0.4 


TVPrT'tir.n riT roc-trvcorH maQQ 

11 CLL. LlVJll VJl 1 V_ __ LUl V_.V_1 llltlji_)i_) 


IMF 


pns rwi rrm 

L,V_i[_. . liJZsl ILJOI 




0.03 


T>1 ot Q 1 "Violri V". ^ V (I'O VI f'1" CI 1 1 All 
111L.LC11 yiv_lv_l |Jv_.l gLlll_,I clLrlvJIl 


i ~\ neon r 

U fJOUI V . 


pn c [77! [70! 


7 

■^prc 


10 -6 


m ot n 1 1 ir'it^/ rHiifi tr. nvp_c>ri r~\ fi~\ mpnt 

llll__U<_llllV_.lLry V_ll__L- L W Ml V_. l__ 111 1L.11111V_.111j 




qprt ISl 

__ C V_- L • 


/Zcj 


5 


tv apt'invi t~iT in nfa lo r^irr'tr' n t/r. nalr^ 

11 CLL, L1U11 Ul I11v_-Iji_1j1_. L.J V_-L. Ll_.tL L VJ IlcLIU 


free 


pnc; _77l l7Kl 






c_ t" v* f" /" \ v* m *_» t" i /"k n o n *H tc-Oiri _-_ f L"' 
3Lal 1UI lllatlUIl cillv.1 ititrtl l_/ctt.rv. 






7th, D 


5 


t n m a 1 otti r h M"^T\r'\T dt \ o f \ \\ ■ i in tnin c\fctf*m 

Lllt.1 llldll L_llll_.lV_.IlV_. V Ul 1L_ V_.U fJclV_.l\. Ill Llllll __ V t) LL. lilt) 


fr 


pn FnTn 

L LJj. IUUJ 


Jth,B 


u.o 


thermal efficiency of feedback in thick systems 


tree 


on 1 ( w 1 


J kin 


o 


l/'ivijrit"!/^ / ^ vi nr in' frnm nr.t unnrlc 

1S_1I1v_.L1v_. V_.IlV_.ltiy 11 Ulll 11U L W lllLlo 




pnc; |(i 1 I If wl 

L.l_|__ . |U -L 1 IU 1 1 




UlJ Kill a 


turbulent velocity of clouds 


free 


pne rvTi r^Ri 


^thr 


D M/^ nr" 2 

U 1V1(V) ^Jv< 


o'qc QiifTt-rT 1 H r>n <ai t ~\r t nfrOQlirnlin Tr.r cf nr fni'iTiatinn 

t_! tViO Lll l<_iv_L_ V_ll_.1101 L y Lilll L_ OllUltl 1 Wl O L<X1 1U1 111<_LL1U11 


y.) ijol_.i v . 


Gprtinn 17 91 

OL_ L_ LlVJll ^^^^^ 


■^hmit 


A/T/-\ i~ii" ^ 

LaJ IV J. (V) U L_ 


r*i*it"i/^_.l fyj.c cnrtarri riir^vjcit'VT" tt.v" riic/^tt 
L.1 1 LrlVjCll Kclo oU.Ilclv_.LJ UtllolLy IUI Ulov_.__ 




pn r^rn 


hot kin. fb 


YES 


switch for heating cold gas by kinetic feedback 


free 


ens 171 1751 






AGN 






-^seed 


1000 M 


seed black hole mass 


theory 


sect. [97T1 


/low J 


0.003 


rate of loss of angular momentum 


free 


eq.[79] 


/jet,0 


1 


efficiency of jet feedback for a 1000 km s -1 halo 


N-bocfy 


eq.[82] 






numerical parameters 






-Mpart 


10 9 M Q 


particle mass 








0.1 Gyr 


numerical interval for the integration 


free 


sect. [22] 



Table 4. Model parameters, with their value adopted in the reference model, brief description, available constraints (independently of 
the model) and reference in the text. Parameters highlighted by a mark are of primary importance. Cosmological parameters are not 
included. 



is fixed by requiring a good match of the power-law part of 
the stellar mass function. 

/th,B, /kin and hot kin. fb: they have only a modest 
influence on the results, so they are left to their reference 
values. 

(To: it influences star formation in small bulges, and then 
the amount of downsizing of the AGN population; it is fixed 
as explained in Fontanot et al. (2006a). 

S t h r : it influences the mass of cold gas in local discs and 
its value is suggested by observations. 

Sii m i t : it influences the history of star formation and 
BH accretion in a rather subtle way; it is used in Fontanot 
et al. (2006a). 

/lowj: it regulates the BH-bulge relation at z = 0, and 
its value depends on whether quasar-triggered winds are 
switched on; see Fontanot et al. (2006a) for details. 

/jct.o: its value is suggested to be unity by N-body sim- 



ulations, and it allows to quench cooling flows and limit the 
mass of elliptical galaxies. 



11 POST PROCESSING 

For each tree that is analyzed, all the variables introduced 
in section 12.31 together with the disc and bulge radii and 
velocities, star formation rates and accretion rate onto black 
holes (at the sampling time, not averaged over the time bin), 
are sampled in intervals of time At and at the final time. 
All these quantities are output at the end, together with 
the merger histories of galaxies, constructed as specified in 
section 13.11 

For each time bin this information is issued for all the 
existing galaxies of a tree. However, this is not sufficient to 
reconstruct the star-formation history of a galaxy compo- 
nent, because disc instabilities, mergers, tidal stripping and 
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destruction events move stars among galaxies and compo- 
nents; as a consequence, the sampled star formation history 
of stars that are formed in a galaxy component does not 
coincide with the star formation history of the stars even- 
tually found in that galaxy component. To reconstruct this 
last quantity, all the events that transfer stellar mass from a 
component or a galaxy to another are recorded and output 
at the end. Star formation histories are then reconstructed 
at the post-processing stage by scrolling the sampled history 
of galaxies and moving stellar mass among the components. 
This is a very quick process, that can be done at the reading 
time. 

The output of MORGANA (star formation histories and 
metallicities of the cold gas) is then given to the spectro- 
photometric code GRASIL (Silva et al. 1998), which is able 
to predict the SEDs of the predicted galaxies from the UV 
to the radio. This way it is possible to construct predictions 
for the luminosity functions at a fixed redshift or in redshift 
intervals, galaxy number counts and galaxy backgrounds. 
The GRASIL code has been tested against many local and 
distant galaxies and is widely used by the community; how- 
ever, it introduces further uncertainties and parameters, so 
we prefer in this paper to discuss only the prediction of MOR- 
GANA before its interfacing with GRASIL. All the details of 
this interfacing will be given in Fontanot et al. (2006b). 

The star formation rates are then reconstructed in time 
bins of width At, which is set to 0.1 Gyr. This sampling is 
fine but for the last bin considered, as the last formed stars, 
that dominate the UV, B and FIR spectra, live for less than 
the sampling bin. To overcome this difficulty, we sample also 
the value of bulge and disc star-formation rates evaluated at 
the end of each time bin. Then, for each galaxy component, 
the time bin corresponding to the time at which the SED is 
computed is split into two parts, of widths 0.09 and 0.01 Gyr. 
To the second part we assign a star-formation rate equal to 
the punctual value at the end of the bin, while to the first 
part we assign a star-formation rate such as the integral over 
the bin gives the total amount of stars. As GRASIL re-samples 
the star-formation histories on a much finer time grid, this 
procedure is accurate as long as the star formation rate does 
not vary strongly on timescales smaller than 10 Myr, which 
is the case for most galaxies. The accuracy of this procedure 
with respect to the use of the complete star-formation rate 
produced by MORGANA during the integration will be shown 
in Fontanot et al. (2006b). 

To compute the AGN activity of galaxies, a sampling 
in a time grid is not sufficient, due to the low duty cycle of 
AGNs. To optimize the statistical sampling of AGNs with- 
out inflating the output files, we save all significant accretion 
events at each integration time-step; this is done whenever 
the punctual value of the accretion rate is larger than a 
minimal value of 1.76 x 10 -3 Mq yr _1 , corresponding to 
a bolometric luminosity of 10 43 erg s _1 . This very detailed 
output adds only ~10 per cent to the total disc space needed 
by a run. For the computation of the luminosity function 
of AGNs, each of these events is treated as an indepen- 
dent event of duration equal to the integration interval of 
time. This procedure is explained in detail in Fontanot et 
al. (2006a). 



12 RESULTS 

As already mentioned in the Introduction, this paper is de- 
voted to a detailed description of the MORGANA code, so in 
this section we present only some of the main results, ob- 
tained with the set of parameters given in table U unless 
otherwise stated. For sake of simplicity we restrict ourselves 
to predictions that do not require the use of a spectropho- 
tometric code. In this way all the uncertainties involved in 
the generation of SEDs for the model galaxies are bypassed. 

All the results are based on a single 512 3 PINOCCHIO run 
of a 150 Mpc box (h = 0.7) with the standard cosmology 
given in the Introduction. The particle mass is 1.0 x 10 9 
Mq; then the smallest considered halo contains 50 particles, 
for a mass of 5.1 x 10 10 M , while the mass of the smallest 
progenitor is 1.0 x 10 10 Mq. As a comparison, the much 
bigger Millennium Simulation (Springel et al. 2005) has a 
slightly higher particle mass, 1.2 x 10 9 Mq for our Hubble 
constant, and its merger histories are reconstructed starting 
from 20 particles, corresponding to a smallest progenitor 
of 2.5 x 10 10 . To test the stability of the results, two more 
boxes have been produced with the same number of particles 
and sizes of 200 and 100 Mpc, for particle masses a factor 
of 2.4 worse (higher) and 3.4 better (lower). The results of 
the stability tests are briefly shown in appendix [B] we will 
highlight in the following which results are more sensitive to 
mass resolution. 

The MORGANA runs have been performed on simple 
PCs. A 150 Mpc run with nearly 2000 trees needs about 
20 hours on a 3GHz PC, the generated output is about 3 
Gb. The bottleneck of the computation is GRASIL, which re- 
quires about five minutes per galaxy, so the computation of 
10000 galaxies in a mock pencil beam survey requires about 
one month of CPU. 

We find that the stellar mass of the typical central 
galaxy contained in the smallest DM halos at z = is 3 x 10 8 
Mq, so the stellar mass function is severely incomplete below 
this limit. 

To limit the size of the output and the computing time, 
the mass function of DM halos is sampled by picking no 
more than 300 halos per mass bin of 0.5 dex (or all halos in 
the bin if they are less than 300). The statistical distribu- 
tions of galaxies are then computed by weighting each galaxy 
by TOtrcc, the inverse of the fraction of the host DM halos 
picked up in the mass bin. For each DM halo all galaxies 
and satellites are computed, so that for a given stellar mass 
the satellites of large DM halos are oversampled with respect 
to the central galaxies, which are numerically dominant at 
all masses. This feature shows up when galaxy properties 
are reported for samples of model galaxies, like in the Tully- 
Fisher relation: at fixed stellar mass, satellite galaxies with a 
low weight Wtroo are more numerous in the plot than central 
galaxies with a high w tr ee- We then further sparse-sample the 
satellites as follows: we compute the average stellar mass of 
central galaxies as a function of the DM halo mass, then as- 
sign to each satellite a weight w ga i equal to the ratio between 
the u; trC c °f the average DM halo hosting a central galaxy 
with the same stellar mass and the Wt rcc of the host DM 
halo. We then select the satellite with a probability equal to 
1/wgai. In this way we obtain a roughly constant number of 
galaxies per logarithmic interval of stellar mass. 

Figure [6] gives, for each selected DM halo evolved to 
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Figure 6. DM halos at z = 0: fraction of baryons in hot halo gas (red triangles, upper left panel), expelled by winds (magenta circles, 
upper right panel), in stars in all galaxies (green open squares, lower left panel), in stars in the central galaxy (blue crosses, lower right 
panel). Averages and rms of the points are shown in all panels. 



2 = 0, the fraction of baryons that are found in hot halo 
gas, in stars of all galaxies (but not in the halo), in stars 
of the central galaxy, and the fraction of baryons ejected by 
winds and not fallen back. These clearly do not exhaust the 
list of baryonic components (halo stars and cold gas compo- 
nents are excluded here), so these fractions do not add up 
to one. However, this figure helps in showing a number of 
important features. There is a clear transition at a DM halo 
mass range of 10 12 — 10 13 Mq: smaller halos lose most of 
their baryons by ejecting them to the IGM and their stellar 
content is dominated by the central galaxy, while larger ha- 
los retain most of their baryons as hot gas and have a small 



fraction of stars (~ 20 per cent) in the central galaxy. In 
both cases the fraction of stars (both total and in the cen- 
tral galaxy) declines quickly with increasing or decreasing 
mass, while a maximum is reached at 10 12 Mq. All these 
features reproduce nicely the trends suggested in the semi- 
nal papers of galaxy formation (see, for instance, White & 
Rees 1978 or Dekel & Silk 1986): small halos are evacuated 
by winds driven by SNe, while in large halos the large cool- 
ing time (aided by AGN feedback) prevents most gas from 
cooling. As a further element, it is worth noting the presence 
of many Milky Way-sized halos with almost no hot gas, in 
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Figure 7. Stellar mass function of galaxies at z = (continuous 
black line), compared with the results of 2dF+2MASS (Cole et 
al. 2001) and SDSS (Bell et al. 2003). The dashed red and dotted 
blue lines give respectively bulge-dominated and disc-dominated 
galaxies, while the dot-dashed magenta line refers to a model 
without quenching of the cooling flows by AGN jets (/j o t,0 = 0). 



Figure 8. Star formation rate density as a function of redshift, 
compared with the data compiled and homogenized by Hopkins 
(2004). Points denote star formation estimates based on rest- 
frame UV (blue squares), [Oil] (green triangles), H« and H/3 
(magenta circles), X-ray, FIR and sub- mm (red stars). 



nice agreement with the lack of X-ray emitting around our 
Galaxy. 

The trends of decreasing efficiency of star formation in 
smaller and larger halos than 10 12 Mq are needed to pro- 
duce the observed low-luminosity slope and high-luminosity 
cutoff of the galaxy LFs. Figure [7] shows the predicted stel- 
lar mass function of galaxies, compared with the data by the 
2dF+2MASS (Cole et al. 2001) and SDSS (Bell et al. 2003) 
surveys. The low-mass slope is nicely reproduced (though 
the result is sensitive to resolution, see appendix [Bj , even 
in the slight steepening below 10 10 Mq . The high-mass cut- 
off is not strong enough, and the biggest ellipticals are too 
massive by at most a factor of two; however, this excess 
may be connected, as proposed by Monaco et al. (2006), to 
the construction of the diffuse stellar component of galaxy 
clusters, and would be fixed by a proper modeling of the 
scattering of stars during mergers (section |4.3[) ; this work is 
in progress. We also report in the figure the resulting mass 
function without AGN feedback, to show that the quenching 
of the cooling flow by AGN jets decreases significantly the 
discrepancy with observations. 

Figure [7] shows also the stellar mass function of bulge- 
and disc-dominated galaxies (these are very similar to the 
mass functions of bulges and discs). In line with what is 
observed, bulges dominate at large masses, while discs are 
more abundant at small masses. However, as noticed also by 
Croton et al. (2006), the mass function of bulges does not 
have a broad peak at ~ 10 10 — 10 11 Mq, as shown by the 
luminosity function of ellipticals; there is a definite excess 
of small bulges. This is likely connected to the excess in 
the predicted number density of 10 10 Mq galaxies at z ~ 1 
reported by Fontana et al. (2006), a problem shared by many 



galaxy formation models. This point will be deepened in 
Fontanot et al. (2006b). 

Figure [5] shows the prediction of the star formation rate 
density as a function of redshift. The data have been col- 
lected and homogenized by Hopkins (2004). The prediction 
is consistent with the data, with some possible underesti- 
mate at z ~ 1. In particular, a very broad peak of star 
formation is predicted to be present at z ~ 3, in agreement 
with the estimates based on sub-mm counts. The level of star 
formation is still high at z ~ 6, only a factor of 2 lower than 
the peak value; however, this prediction depends sensitively 
on the mass resolution. Most of this star formation takes 
place in discs, triggered by the strong cooling flows at high 
redshift, while mergers dominate the strongest starbursts; 
this will be deepened in Fontanot et al. (2006b). 

Figure||]shows the cosmic stellar mass density as a func- 
tion of redshift. Though this quantity is related to the cosmic 
star formation rate (it is simply its integral in time), it is 
compared with a completely different set of data (Brinch- 
mann & Ellis 2000; Cole et al. 2001; Cohen 2002; Dickinson 
et al. 2003; Fontana et al. 2003, 2004; Drory et al. 2004; 
Gwinn & Hartwick 2005) . The data in this figure are not ho- 
mogenized for the different mass limits, so each point should 
be considered as a lower limit. The agreement with data is 
again good from z ~ 4 to 0. A more refined analysis of the 
build-up of the stellar mass is reported in Fontana et al. 
(2006) and Fontanot et al. (2006b). 

Figure [10] (left panel) shows the baryonic Tully-Fisher 
relation for the disc-dominated galaxies, where disc veloci- 
ties Vd are computed at the optical radius (3.27?d) and Mu 
is the total disc mass. This relation is compared with that 
obtained from the universal rotation curve (Persic, Salucci & 
Stel 1996, updated by Yegorova & Salucci 2006) ; the intrin- 
sic scatter around this relation is reported by the authors 
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Figure 10. Baryonic Tully-Fisher relation for disc-dominated galaxies at z = 0, compared with the relation obtained from the Universal 
Rotation Curve (Persic, Salucci & Stel 1996; Yegorova & Salucci 2006), shown as a shaded area. Dots and circles correspond respectively 
to discs with a gas fraction smaller and higher than 1 per cent. Continuous lines give the average (thick line) and average ±rms (thin 
lines) of the circles. Left panel: standard model; right panel: model with concentration scaled to c n f w = 4 for 10 12 Mq DM halos at 

2 = 0. 
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Figure 9. Stellar mass density as a function of redshift, compared 
with the data by Brinchmann & Ellis (2000), Cole et al. (2001), 
Cohen (2002), Dickinson et al. (2003), Fontana et al. (2003; 2004), 
Drory et al. (2004), Gwinn & Hartwick (2005). The data are not 
homogenized for the different sample completeness. 



to be ~0.2 dex. As a first point, there is a tail of rapidly 
spinning discs that are not seen in the data. These discs are 
formed at high redshift (when DM halos were denser and 
discs more compact), then become satellites so that no fur- 
ther infall on the disc takes place. Consequently, they are 
compact and gas-poor, and thus can hardly be classified as 



spirals. This is demonstrated in the figure, where all objects 
with a gas fraction lower than 1 per cent are denoted as 
small dots; most if not all outliers are then removed by ap- 
plying this selection. As a second point, the gas-rich discs 
follow a Tully-Fisher relation parallel to the observed one 
but with a higher velocity by 0.1 dex, or 25 per cent, in 
Vd. (Notably, galaxies follow the same Tully-Fisher relation 
independently of their surface brightness.) We have verified 
that no realistic combination of feedback parameters leads 
to a better zero-point of the relation. This disagreement is 
more likely related to the shape of the halo; for instance, 
Mo & Mao (2000) suggest that a concentration as low as 
c n fw = 4 is required to fit the Tully-Fisher relation. To test 
this hypothesis we then scaled all concentrations by a con- 
stant factor so as to take a value of 4 for a 10 12 M© halo at 
z = 0. The resulting relation (figure flOl right panel) shows a 
much better agreement with the observed one (this change 
does not affect much the other galaxy properties, but el- 
lipticals result larger and less dense as well). We then con- 
clude that the disagreement in the zero-point of the baryonic 
Tully-Fisher relation is likely due to the inner profile of the 
DM halo more than to feedback. As a third and final point, 
the scatter around the Tully-Fisher relation is found to be 
~0.25 dex, possibly higher than what is suggested by data. 

ft must be stressed that this model was run without 
the computation of adiabatic contraction of the halo (sec- 
tion [671} • We have verified that this further process, besides 
slowing down the computation significantly, makes the discs 
slightly more compact, without influencing drastically the 
other results; in fact, taking into account the presence of 
a bulge leads to a significant contraction in many discs, so 
that the adiabatic contraction of the halo does add much to 
this effect. 

Figure [11] shows the surface densities of cold gas versus 
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Figure 11. Schmidt-Kennicutt law for the disc (blue triangles) 
and bulge (red circles) components. The shaded area gives the 
observed relation reported by Kennicutt (1998). 

star formation rate for discs and bulges, compared with the 
Schmidt-Kennicutt law (Kennicutt 1998). By construction, 
bulges stay exactly on the average relation, while discs fol- 
low the star-formation timescale of equation 1581 which, as 
anticipated, gives a relation compatible with the observed 
one in virtue of the correlation between cold gas fraction 
and surface density. The agreement is very good both in 
terms of zero-point (apart from a marginal underestimate 
that could be easily fixed by a better tuning) and in terms 
of scatter. 

Figure [TZ] shows the prediction of the mass function of 
HI at 2 = 0, compared with that obtained from the HIPASS 
sample (Zwaan et al. 2003). The statistical error of the ob- 
servational mass function (well fit by a Schechter function) 
is so small for this sample that we do not even report it 
in the figure; conversely, the transformation from HI to to- 
tal cold gas mass is rather uncertain. Following Fukugita & 
Peebles (2004), the HI and molecular gas densities at z = 
amount to Q, H i = (3.5 ±0.8) x 10" 4 (Zwaan et al. 2003) and 
Qh 2 = (1.6 ± 0.4) x 10" 4 (Keres, Yun & Young 2003), so 
that, taking into account the abundance of helium (another 
1.38 factor) the fraction of cold gas results a factor of two 
higher than that of HI; this is the conversion factor used 
in figure [12] With the standard choice of parameters MOR- 
GANA does not reproduce this mass function either at small 
or at large masses. In appendix [B] we show that the cutoff 
of this function is remarkably sensitive to mass resolution 
and decreases with decreasing particle mass, so that at this 
stage we decide not to give credit to the discrepancy at large 
masses. At smaller masses, we find that the amount of gas is 
sensitive both to the threshold for star formation E tnr intro- 
duced in section FT21 and to the assumed star-formation law. 
Given the good reproduction of the Schmidt-Kennicutt law 
(even with a marginal underestimate of the star-formation 
rate, which would imply an increase of gas masses) , the dis- 
crepancy in the amount of cold gas cannot be due to an error 
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Figure 12. Mass function of HI gas (black continuous line), com- 
pared with the analytic fit of the results of the HIPASS sample 
(Zwaan et al. 2003; green thick line), assuming that the correc- 
tion for helium and molecular gas amounts to a factor of 2. The 
statistical error of the observed mass function is very small. The 
dashed red and dot-dashed blue lines give the results of models 
with a threshold for star formation (E tnr = 10 Mq pc -2 ); in the 
blue dot-dashed model DM halo concentrations c n f w are scaled 
so as to assume a value of 4 for a 10 12 Mq halo at z = 0. 



in the star- formation time-scale at a given £ C oid,D, but can 
be due to an error in the size of discs, which influences the 
value of E co id.D and then of i*,D- Such an error is also sug- 
gested by the discrepancy, reported above, in the zero-point 
of the baryonic Tully-Fisher. We then show two more mod- 
els in figure WZl both with E t hr = 5 Mq pc -2 , the second 
one with concentrations scaled to take a value c n f w = 4 for a 
10 12 DM halo at z = 0. Clearly, adding a threshold for star 
formation helps in increasing the amount of cold gas (un- 
der the assumption that the correction for the HI gas given 
above applies also beyond the threshold for star formation) 
but cannot solve the discrepancy, while the c n f w = 4 case 
leads even to a satisfactory prediction of the low-end of the 
gas mass function. We conclude that the mass function of 
cold gas gives a very important fine constraint to the model, 
and that it gives further support to the idea that the NFW 
profile of DM halos leads to too compact discs. However, the 
uncertainty in the relation between HI and cold gas must be 
better taken into account before reaching strong conclusions. 

Figure [13] shows the evolution of the cold gas density, 
compared with the HIPASS point (shown with and without 
the correction for molecular gas and helium) and the data 
from DLA systems (Peroux et al. 2005; Rao et al. 2005; 
Prochaska, Herbert-Fort & Wolfe 2005), to which no correc- 
tion has been applied. The same three models as in figure [TZI 
are shown. All the models show the same trend of a broad 
peak at z ~ 2 — 4 and a decrease by a factor of ~ 5 at z < 2; 
the Cntw = 4 model gives a remarkably high normalization. 
The data themselves are characterized by large errorbars 
and by a discrepancy between the Prochaska et al (2005) 
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Figure 13. Evolution of the cold gas density, compared with 
the HIPASS point at z = (with and without the correction for 
helium and molecular gas) and to estimates of HI gas density 
from DLAs, obtained by Peroux et al. (2005), Rao et al. (2005) 
and Prochaska, et al. (2005). Model lines are as in figure 1121 
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Figure 15. Stellar mass versus stellar surface density for disc- 
dominated (blue triangles) and bulge-dominated (red circles) 
galaxies in the model. The thick and thin continuous lines give 
the average and average ± rms of the galaxies. The shaded area 
reports the observational result of Kauffmann et al. (2003). 




Figure 14. Gas content of disc-dominated (blue triangles) and 
bulge-dominated (red triangles) galaxies as a function of stellar 
mass. 



and Peroux et al. (2005) datapoints at z ~ 2.5. No strong 
conclusion can then be drawn from this figure, apart from a 
rough qualitative agreement. 

FigurelTilshows. for the standard model, how cold gas is 
distributed in disc-dominated or bulge-dominated galaxies. 
Regarding the latters, while small ellipticals have a low gas 
content, the most massive objects have a significant gas load; 
at variance with what appears in the plot, they do not dom- 



inate the high end of gas mass function because they are 
hosted with the rare but well sampled massive DM halos. 
This excess of cold gas does not contribute much in terms 
of mass, but makes these galaxies too blue, thus highlight- 
ing that our quenching of the cooling flow by AGN feedback 
is not strong enough. This point will be addressed in more 
detail in a forthcoming paper. 

Figure [151 shows the galaxies in the stellar mass-surface 
density plane. Spiral discs and bulges tend to occupy dif- 
ferent regions of the plane; the scatter is however so strong 
that, as observed, no clear bimodality emerges from the dis- 
tribution of these galaxies. The region occupied by SDSS 
galaxies in this plot (Kauffmann et al. 2003) is highlighted 
in the figure. A decrease of the surface brightness at masses 
< 10 10 Mq is obtained, in good agreement with the SDSS 
data; however, model galaxies are much more scattered. 

Figure \W\ shows the bulge-dominated galaxies in the 
stellar mass-i?B plane, compared with the observed relation 
obtained from the data of Marconi & Hunt (2003), Rb = 
4.9(M B /10 n M ) O<35 kpc (see also Monaco & Fontanot 
2005), with a scatter of 0.3 dex. The data lie nicely within 
the observed range, with some flattening at small masses 
which is consistent with the findings of Graham et al. (2006) . 
We do not show a prediction of the fundamental plane at this 
stage, because the virial theorem is implicit in our relation 
between mass, radius and velocity dispersion, and the com- 
putation of M/C ratios need a spectro-photometric code to 
be computed. 

Figure [T7] shows the predicted metallicity of bulges, 
discs and hot halo gas. Elliptical galaxies show a mass- 
metallicity relation that steepens at masses smaller than 
ll 10 Mq, spirals show a similar though weaker trend. These 
metallicities are compared with the results of Nelan et al. 
(2005), relative to the cluster elliptical galaxies of the NOAO 
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fundamental plane survey, and Gallazzi et al. (2005), rela- 
tive to SDSS galaxies. Both datasets suggest an increase of 
metallicity with stellar mass, which confirms our qualita- 
tive prediction, but they differ remarkably in the normaliza- 
tion. This reflects the intrinsic difficulties in estimating ages 
and metallicities in large samples of galaxies. Taking this 
into account, we notice that the region of masses lower than 
~ 3 x 10 10 Mq and metallicities lower than 0.005 is popu- 
lated by SDSS galaxies (presumably by disc-dominated ob- 
jects) but not by our model. Taking this result at face value, 
these data suggest that there may be a need of feedback in 
small-mass discs. In Figure [17] we also show the predicted 
metallicity of hot halo gas, which is enriched to a level of 
roughly Zq/3 in clusters, which raises in galaxy groups and 
then is dominated by scatter at smaller DM halo masses 
(which are also almost devoid of hot gas). This is in quali- 
tative agreement with the observed trend (Baumgartner et 
al. 2005) of an increase of metallicity from galaxy clusters 
to groups of 2-3 keV (~ 10 14 M©), possibly followed by a 
decrease at smaller masses. 



Figure 16. Half-mass radii of elliptical galaxies as a function of 
stellar mass. Points with error bars give the average and rms in 
bins of the model galaxies. The lines give the the average and rms 
values of the data obtained from Marconi & Hunt (2003). 
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Figure 17. Left panel: metallicity of bulge-dominated (red 
points) and disc-dominated (blue triangles) galaxies, compared 
to the average and rms range of data of SDSS galaxies (Gallazzi 
et al. 2005, black dashed lines) and of the NOAO survey clus- 
ter ellipticals (Nelan et al. 2005, green dotted lines). Right panel: 
predicted metallicity of the hot halo component; in this panel, the 
continuous thick and thin cyan lines give the average and rms of 
the points. 



13 DISCUSSION AND CONCLUSIONS 

We have presented and described in detail the code MOR- 
GANA for the formation and evolution of galaxies and AGNs. 
The most relevant features of MORGANA have been men- 
tioned in the Introduction: it attempts, through its original 
modeling of cooling, star formation, feedback, galactic winds 
and super-winds, AGN activity and AGN feedback, to move 
from a phenomenological description of galaxy formation, 
based on simple scalings with the properties of the host DM 
halo, toward a fully physically motivated one. The numeri- 
cal integration of the mass, energy and metal flows allows a 
wide set of physical processes to be straighforwardly imple- 
mented, and the multi-phase description of the inter-stellar 
and intra-cluster media, although used only in a limited way, 
is a step forward towards a more realistic description of the 
complex physics involved. 

Despite many significant technical differences, the pre- 
dictions of MORGANA at z = are in line with that of the 
other semi-analytical or N-body models: most of the figures 
shown in section [121 will be no real surprise for most experts 
of galaxy formation. Besides, a qualitative comparison of our 
results with the N-body ones of Tornatore et al. (2003) has 
shown in many cases the same trends, as for instance the for- 
mation of too massive ellipticals in clusters or the inability 
of stripping to produce a sufficient number of halo stars (see 
also Monaco et al. 2006). Similarly, the hidden dependence 
of the results on mass resolution and the lack of a proper 
convergence is shared by the two methods. This implies that 
the field is reaching an interesting level of maturity, so that 
an increase in the level of sophistication in these models is 
justified. 

In the Introduction we stated that we do not regard this 
model as a "theory of everything" for galaxies, but simply as 
a powerful tool to understand the complex nature of galax- 
ies and to bridge in a realistic way the physical processes 
in play with the observations that can constrain them. It is 
then important to focus on the discrepancies with observa- 
tions found in this paper, and to the insights on the physical 
processes that they provide. 
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The exponential cutoff of the luminosity function is only 
roughly reproduced, which implies that the quenching of the 
cooling flow by the AGN is not well modeled; a stronger in- 
dication is given by the excess of cold gas in large ellipticals. 
This is no surprise given the poor level of understanding of 
this process. Our quenching is performed by allowing gas to 
cool, form stars and accrete onto the BH, at variance with 
the choice of switching cooling off whenever some fiducial cri- 
terion is satisfied, like e.g. in Hatton et al (2003), Bower et 
al. (2006) or Croton et al. (2006), a choice that leads to bet- 
ter fits to the luminosity functions and correct red colours 
for the most massive galaxies. Besides, if the build-up of 
massive galaxies by consecutive mergers is slowed down by 
the creation of the diffuse stellar component in galaxy clus- 
ters, as proposed by Monaco et al. (2006), then the high end 
of the stellar mass function (but not galaxy colours) may be 
determined by gravitational processes more than by AGN 
feedback. This suggests that the longly known cutoff of the 
stellar mass function may still give surprises in the future. 

Another well-known discrepancy arises in reproducing 
the zero-point of the baryonic Tully-Fisher relation. Our re- 
sults support strongly the idea that this problem is not due 
to an unsuitable modeling of feedback but to the shape of 
the inner density profile of DM halos; the cusped NFW halos 
produce discs that are too compact (with Vb values higher 
by 25 per cent), while by scaling the halo concentration so 
as to assume the low value c n f w = 4 for a 10 12 Mq DM halo 
at z = 0, as suggested by Mo & Mao (2000), it is possible to 
recover the zero-point of the Tully-Fisher relation. Besides, 
the formation of galaxy discs is based on the (reasonable but 
still not fully demonstrated) assumptions of conservation of 
angular momentum and formation of exponential discs, and 
our model neglects important problems related to the distri- 
bution of angular momentum within the DM halos. Given 
the complexity of this topic, it is remarkable that the ob- 
served Tully-Fisher relation is wrong by only 0.1 dex under 
the simplest assumptions. 

The discrepancy in disc sizes is reflected also in the 
mass function of cold gas; model disc galaxies have correct 
stellar masses and stay on the correct Schmidt-Kennicutt 
law, but show a deficit of cold gas mass. This lack is partially 
fixed by inserting a threshold for star formation, but only 
the formation of less compact discs (hosted in flatter DM 
halos) can solve the discrepancy. Our analysis also shows the 
potentiality of using the mass function of cold gas in galaxies 
to test galaxy formation models, although the uncertainty in 
the relation between HI and total cold gas hampers strong 
conclusions. 

Another point of disagreement with data is given by the 
excess of small bulges at z — 0. This is related to the excess 
of 10 10 Mq galaxies at z > 1 (Fontana et al. 2006), and to 
the intrinsic difficulty of "downsizing" the AGN population 
(Fontanot et al. 2006a). The lack of low-metallicity discs 
may add to this evidence. These connected problems point to 
some missing feedback sources that acts in disc-dominated 
galaxies at high redshift (those that later merge to form 
bulges) . 

All these discrepancies point to the need of a deeper 
understanding of the complex process of galaxy formation, 
that will necessarily need a focus on the details of the various 
processes involved. In this regard, observations of larger and 
larger samples of galaxies will need to be complemented by 



detailed observations of objects where feedback is at work. 
To bridge together the two sets of evidences it will be neces- 
sary to use galaxy formation models that, on the one hand, 
are able to generate predictions for large samples of ob- 
jects and, on the other hand, contain a sophisticated enough 
treatment of gas dynamics to make predictions on the de- 
tailed properties of galaxies. MORGANA is a contribution in 
this direction. 



ACKNOWLEDGMENTS 

We thank Lucia Ballo, Stefano Borgani, Gabriele Ces- 
cutti, Gabriella De Lucia, Carlos Frenk, Gianluigi Granato, 
Giuseppe Murante, Ezio Pignatelli, Paolo Salucci, Laura 
Silva and Tom Theuns for many discussions, and Stefano 
Cristiani for his encouragement. The anonymous referee has 
been very helpful in improving the presentation of the pa- 
per. PM thanks the Institute for Computational Cosmology 
of Durham for hospitality. 



REFERENCES 

Arnaboldi, M., Freeman, K.C., Okamura, S., et al., 2003, 
AJ, 125,514 

Balland, C, Devriendt, J.E.G., Silk, J., 2003, MNRAS, 
343, 107 

Baugh, CM., Lacey, C.G., Frenk, C.S., Granato, G.L., 
Silva, L., Bressan, A., Benson, A.J., Cole, S., 2005, MN- 
RAS, 356, 1191 

Baumgartner, W.H., Loewenstein, M., Horner, D.J., 
Mushotzky, R.F., 2005, ApJ, 620, 680 
Bell, E.F., Mcintosh, D.H., Katz, N., Weinberg, M.D., 
2003, ApJS, 149, 289 

Benson, A.J., Bower, R.G., Frenk, C.S., Lacey, C.G., 
Baugh, CM., Cole, S., 2003, ApJ, 599, 38 
Benson, A. J., Lacey, C.G., Baugh, CM., Cole, S., Frenk, 
C.S., 2002, MNRAS, 333, 156 

Bianchi, S., Cristiani, S., Kim, T.-S., 2001, A&A, 376, 1 
Bond J.R., Cole S., Efstathiou C, Kaiser N., 1991, ApJ, 
379, 440 

Bower, R.G., Benson, A. J., Malbon, R., Helly, J.C., Frenk, 

C.S., Baugh, CM., Cole, S., La cey, C.G., 2006, submitted 

to MNRAS ( |astro-ph/0511338| ) 

Brinchmann, J., Ellis, R.S., 2000, ApJL, 536, L77 

Bullock, J.S., Dekel, A., Kolatt, T.S., Kravtsov, A.V., 

Klypin, A.A., Porciani, C, Primack, J.R., 2001, ApJ, 555, 

240 

Buchert T., Ehlers J., 1993, MNRAS, 264, 375 
Cattaneo, A., Dekel, A., Devriendt, J., Guiderdoni, B., 
Blaizot, J., 2006, submitted to MNRAS ( |astro-ph/0601295[ ) 
Cavaliere, A., Vittorini, V., 2000, ApJ, 543, 599 
Christodoulou, D.M., Shlosmann, I., Tohline J.E., 1995, 
ApJ, 443, 551 

Cimatti, A.; Pozzetti, L.; Mignoli, M., et al., 2002, A&A, 
391, LI 

Cohen, J., 2002, ApJ, 567, 672 

Cole, S., Lacey, C.G., Baugh, CM., Frenk, C.S., 2000, MN- 
RAS, 319, 168 

Cole, S., Norberg, P., Baugh, CM., et al., 2001, MNRAS, 
326, 255 



© RAS, MNRAS OOO.ITU351 



32 Monaco, Fontanot & Taffoni 



Cox, A.N., 2000, Allen's astrophysical quantities, 4th ed. 
(New York: AIP Press). Pag. 572 

Croton, D.J., Springel, V., White, S.D.M., et al., 2006, MN- 
RAS, 365, 11 

Dekel, A., Birnboim, Y., 2006, MNRAS, 368, 2 

Dekel, A., Silk, J., 1986, ApJ, 303, 39 

De Grandi, S., Molendi, S., 2002, ApJ, 567, 163 

De Lucia, G., Springel, V., White, S.D.M., Croton, D., 

Kauffmann, G., 2006, MNRAS, 366, 499 

Dickinson, M., Papovich, C, Ferguson, H.C., Budavari, T., 

2003, ApJ, 587, 25 

Dib, S., Bell, E., Burkert, A., 2006, ApJ, 638, 797 

Di Matteo, T., Croft, R.A.C., Springel, V., Hernquist, L., 

2003, ApJ, 593, 56 

D'Onghia, E., Burkert, A., Murante, G., Khochfar, S., 
2006, submitted to MNRAS < |astro-ph/0602005[ ) 
Drory, N., Bender, R., Feulner, G., Hopp, U., Maraston, 
C, Snigula, J., Hill, G.J., 2004, ApJ, 608, 742 
Efstathiou, G., Lake, G., Negroponte, J., 1982, MNRAS, 
199, 1069 

Eisenstein, D.J., Zehavi, I., Hogg, D.W., et al., 2005, ApJ, 
633, 560 

Eke, V.R., Navarro, J.F., Steinmetz, M., 2001, ApJ, 554, 
114 

Feldmeier, J. J., Ciardullo, R., Jacoby, G.H., Durrell, P.R., 
2003, ApJS, 145, 65 

Fontana, A., Donnarumma, I., Vanzella, E., et al., 2003, 
ApJL, 594, L9 

Fontana, A., Pozzetti, L., Donnarumma, I., et al., 2004, 
A&A, 424, 23 

Fontana, A., Salimbeni, S.,Grazian, A., et al., 2006, A&A, 
in press ( |astro- ph/0609068| 

Fontanot, F., Monaco, P., Cristiani, S., Tozzi, P., 2006a, 
MNRAS, in press ( |astro-ph/0 609823 ) 

Fontanot, F., Monaco, P., Silva, L., et al., 2006b, in prepa- 
ration 

Genzel, R., Tacconi, L., Eisenhauer, F, et al., 2006, Nature, 
442, 786 

Fraternali, F., Cappi, M., Sancisi, R., Oosterloo, T., 2002, 
ApJ, 578, 109 

Fukugita, M., Peebles, P.J.E., 2004, ApJ, 616, 643 
Gallazzi, A., Chariot, S., Brinchmann, J., White, S.D.M., 
Tremonti, C.A., 2005, MNRAS, 362, 41 
Gal- Yam, A., Maoz, D., Guhathakurta, P., Filippenko, 
A.V., 2003, AJ, 125, 1087 

Governato, F., Mayer, L, Wadsley, J., Gardner J. P., Will- 
man, B., Hayashi, E., Quinn, T., Stadel, J., Lake, G., 2004, 
ApJ, 607, 688 

Governato, F., Willman, B., Mayer, L., Brooks, A., Stinson, 
Valenzuela, O., Wadsley, J., Quinn, T., 2006, submitted to 



MNRAS (astro-ph/0602351) 



Graham, A.W., Merritt, D., Moore, B., Diemand, J., 
Terzic, B., 2006, AJ, in press ( astro-ph/ 0608614) 
Granato G.L., De Zotti G., Silva L., Bressan A., Danese 
L., 2004, ApJ, 600, 580 

Granato G.L., Silva L., Monaco P., Panuzzo P., Salucci P., 
De Zotti G., Danese L., 2001, MNRAS, 324, 757 
Gwyn, S.D.J. , Hartwick, F.D.A., 2005, AJ, 130, 1337 
Hatton S., Devriendt J. E.G., Ninin S., Bouchet F.R., 
Guiderdoni B., Vibert D., 2003, MNRAS, 343, 75 
Hopkins, A.M., 2004, ApJ, 615, 209 

Rang, X., Jing, Y.P, Mo, H.J., Borner, G., 2005, ApJ, 631 



21 

Kauffmann, G., Colberg, J.M., Diaferio, A., White, S.D.M., 
1999, MNRAS, 303, 188 

Kauffmann, G., Heckman, T.M., White, S.D.M., et al., 
2003, MNRAS, 341, 54 

Kazantzidis, S., Mayer, L., Colpi, M., Madau, P., Debat- 
tista, V.P., Wadsley, J., Stadel, J., Quinn, T., Moore, B., 
ApJ, 623, L67 

Kennicutt, R.C.Jr., 1989, ApJ, 344, 685 
Kennicutt, R.C.Jr., 1998, ApJ, 498, 541 
Keres, D., Katz, N., Weinberg, D.H., Dave, R., 2005, MN- 
RAS, 363, 2 

Keres, D., Yun, M.S., Young, J.S., 2003, ApJ, 582, 659 
Knop, R.A., Aldering, C, Amanullah, R., et al., 2003, ApJ, 
598, 102 

Kravtsov, A.V., 2003, ApJ, 590, LI 

Lacey C, Cole S., 1993, MNRAS, 262, 627 

Li, Yun, Mo, H.J., van den Bosch, F.C., 2005, submitted 

to MNRAS ( |astro-ph/0510372[ ) 

Lia, C, Portinari, L., Carraro, G., 2002, MNRAS, 330, 821 
Mac Low, M.-M., 2003, in "Simulations of magnetohydro- 
dynamic turbulence in astrophysics", eds. E. Falgarone & 
T. Passot, Springer, pag. 182 

Mailer, A.H., Bullock, J.S., 2004, MNRAS, 355, 694 

Marconi A., Hunt L.K., 2003, ApJ, 589, 21 

Mathis, H., Lemson, G., Springel, V., Kauffmann, G., 

White, S.D.M., Eldar, A., Dekel, A., 2002, MNRAS, 333, 

739 

Menci, N., Cavaliere, A., Fontana, A., Giallongo, E., Poli, 
F., Vittorini, V., 2004, ApJ, 604, 12 

Merloni, A., Heinz, S., Di Matteo, T., 2003, MNRAS, 345, 
1057 

McNamara, B.R., Nulsen, P.E.J., Wise, M.W., Rafferty, 
D.A., Carilli, C, Sarazin, C.L., Blanton, E.L., 2005, Na- 
ture, 433, 45 

Mo, H.J., Mao, S., White, S.D.M., 1998, MNRAS, 295, 319 

Mo, H.J., Mao, S., 2000, MNRAS, 318, 163 

Monaco, P., 2004a, MNRAS, 352, 181 (M04) 

Monaco, P., 2004b, MNRAS, 354, 151 

Monaco, P., Fontanot, F., 2005, MNRAS, 359, 283 

Monaco, P., Murante, G., Borgani, S., Fontanot, F., 2006, 

ApJ Letters, in press (astro-ph/0610045 ) 

Monaco P., Salucci P., Danese L., 2000, MNRAS, 311, 279 

Monaco, P., Theuns, T., Taffoni, G., Governato, F., Quinn, 

T., Stadel, J., 2002a, ApJ, 564, 8 

Monaco P., Theuns T., Taffoni G., 2002b, MNRAS, 331, 
587 

Moore, B., Katz, N., Lake, G., Dressier, A., Oemler, A. Jr., 
1996, Nature, 379, 613 

Moutarde F.,Alimi J. M. , Bouchet F. R.,Pellat R., Ramani 
A., 1991, ApJ, 382, 377 

Murante, G., Arnaboldi, M., Gerhard, O., et al., 2004, ApJ, 
607, L83 

Nagashima, M., Lacey, C.G., Okamoto, T., Baugh, CM., 
Frenk, C.S.; Cole, S., 2005, MNRAS, 363, L31s 
Nagashima, M., Yahagi, H., Enoki, M., Yoshii, Y., Gouda, 
N., 2005, ApJ, 634, 26 

Navarro J.F., Frenk C.S., White S.D.M., 1995, MNRAS, 
275, 720 

Nelan, J.E., Smith, R.J., Hudson, M.J., Wegner, G.A., 
Lucey, J.R., Moore, S.A.W., Quinney, S.J., Suntzeff, N.B., 
ApJ, 2005, 632 137 



© RAS, MNRAS OOO.ITH351 



MORGANA 33 



Ostriker, J. P., Bode, P., Babul, A., 2005, ApJ, 634 964 
Pearce, F. R., Jenkins, A., Frenk, C.S., White, S.D.M., 
Thomas, P.A., Couchman, H.M.P., Peacock, J.A., Efs- 
tathiou, G., 2001, MNRAS, 326, 649 

Percival W.J., Sutherland W.J., Peacock J.A., et al., 2002, 
MNRAS, 337, 1068 

Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, 
ApJ, 517, 565 

P/'eroux, C, Dessauges-Zavadsky, M., D'Odorico, S., Kim, 
T.S., McMahon, R.G., 2005, MNRAS, 363, 479 
Persic, M., Salucci, P., Stel, P., 1996, MNRAS, 281, 27 
Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, 
B.P, 1992, Numerical recipes in FORTRAN, Cambridge: 
University Press 

Prochaska, J., Herbert-Fort, S., Wolfe, A., 2005, ApJ, 635, 
123 

Rao, S.M., Prochaska, J.X., Howk, J.C, Wolfe, A.M., 2005, 
AJ, 129, 9 

Recchi, S., Matteucci, F., D'Ercole, A., Tosi, M., 2002, 
A&A, 384, 799 

Refregier, A., 2003, ARA&A, 41, 645 

Rosati, P., Borgani, S., Norman, C, 2002, ARA&A, 40, 539 

Schaye, J., 2004, ApJ, 609, 667 

Schmidt M., 1959, ApJ, 129, 243 

Silva, L., Granato, G.L., Bressan, A., Danese, L., 1998, 
ApJ, 509, 103 

Somerville, R., Primack, J.R., 1999, MNRAS, 310, 1087 
Spergel, D.N., Bean, R., Dore, O., et al., 2006, submitted 
to ApJ ( |astro-ph/0603449[ ) 

Spergel, D.N., Verde, L., Peiris, H. V., et al., 2003, ApJS, 
148, 175 

Springel, V., Hernquist, 2003, MNRAS 339, 289 
Springel, V., White, S.D.M., Jenkins, A., et al., 2005, Na- 
ture, 435, 629 

Steinmetz, M., Navarro, J.F., 2002, New A, 7, 155 
Sutherland R., Dopita M., 1993, ApJS, 88, 253 
Suto, Y., Sasaki, S., Makino, N., 1998, ApJ, 509, 544 
Taffoni G., Mayer L., Colpi M., Governato F., 2003, MN- 
RAS, 341, 434 

Taffoni G., Monaco P., Theuns T., 2002, MNRAS, 333, 623 
Thomas D., Maraston, C, Bender, R., de Oliveira, CM., 
2005, ApJ, 621, 673 

Toft, S., Rasmussen, J., Sommer-Larsen, J., Pedersen, K., 

2002, MNRAS, 335, 799 

Tormen, G., 1997, MNRAS, 290, 411 

Tornatore, L., Borgani, S., Springel, V., Matteucci, F., 
Menci, N., Murante, G., 2003, MNRAS 342, 1025 
Umemura M., 2001, ApJ, 560, L29 

Viel, M., Haehnelt, M.G., Springel, V., 2004, MNRAS, 354, 
684 

Voit, G.M., Donahue, M., 2005, ApJ, 634, 955 
Volonteri, M., Haardt, F., Madau, P., 2003, ApJ, 582, 559 
Warren, M.S., Quinn, P.J., Salmon, J.K., Zurek, W.H., 
1992, ApJ, 399, 405 

Weinberg D.H., Hernquist L., Katz N., 2002, ApJ, 571, 15 
White, S.D.M., Frenk, C.S., 1991, ApJ, 379, 52 
White, S.D.M., Rees, M.J., 1978, MNRAS, 183, 341 
Wu, K.K.S., Fabian, A.C., Nulsen, P.E.J., 2000, MNRAS, 
318, 889 

Yegorova, Y., Salucci, P., 2006, in preparation 
Young, P.J., 1976, AJ, 81, 807 



Zhao, D.H., Jing, Y.P., Mo, H.J., Borner, G., 2003, ApJ, 
597, L9 

Zibetti, S., White, S.D.M., Schneider, D.P., Brinkmann, J., 
2005, MNRAS, 358, 949 

Zwaan, M.A., Staveley-Smith, L., Koribalski, B.S., et al., 
2003, AJ, 125, 2842 



APPENDIX A: MERGING AND 
DESTRUCTION TIMES FOR GALAXIES 

The merging and destruction times are computed using the 
physically motivated analytic fits (accurate to ~ 15 per cent) 
of Taffoni et al. (2003) to the merging and destruction times 
found using N-body simulations and analytical models. We 
use here a version of the analytic fits that differs slightly 
from the one proposed in that paper in the values of the 
fitting parameters: 

Let the merging halos have mass, circular velocities, 
virial radii and concentrations respectively Mh, Vh, rn, 
ch (main halo) and Ms, Vs, rs, cs (satellite). Let / sat = 
Ms,o/Mh be the ratio between the initial mass of the satel- 
lite and the halo mass (we recall that the halo mass includes 
the satellite), and let e and x c be the orbital parameters of 
the satellite (see section |4~T|| . In the case of a rigid satellite 
that suffers no mass loss, the dynamical friction time (the 
time required by an orbit to decay) is: 

Trigid = 0.46 



GMs,o 



(1.7265 + 0.0416c H ) 



ln(l + l// sa t) 



(Al) 



while for a live satellite, i.e. a satellite that suffers significant 
mass loss by tidal stripping, the dynamical friction time is: 



THve = T^rrz X 



GMs fi 

x [0.25 (ch/cs) 6 - 0.07(c s /c H ) + 1.123] 
x [0.4 + (e - 0.2) Q] 



(A2) 



where 



B = -0.0504 + 0.3355x c + 0.3281a; 

2 



(A3) 

C = 2.151 - 14.176x c + 27.383^ (A4) 
Q = 0.9 + 10 8 (/ sat - 0.0077/(1 - 1.08a; c ) - 0.0362) 6 

(12.84 + 3.04z c - 23.4x 2 c ) (A5) 

The general expression for the dynamical friction time is 
obtained by interpolating between the two expressions: for 
/sat > 0.1: 

Tdf = frigid ; (A6) 
for 0.08 < < 0.1: 

/sat - 0.08 



Tdf = T, 



igid I " 



+ r livc[ l-h^) ; (A7 ) 



"live ^1 



0.02 



(A8) 



0.02 

for 0.007 < /sat ^ 0.0 

Tdf = Tli vc ; 

for / sat s? 0.007: 

r df = max (2.1r rigide °- 475 ( 1 - tanh ( 10 - 3/ - 3 -^^c)) 5Tiiv ^ (A9) 

For the destruction times we follow closely Appendix B 
of Taffoni et al. (2003). 
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APPENDIX B: NUMERICAL STABILITY 

In this Appendix we show how some results are sensitive 
to the mass resolution of the merger trees. This is not to 
be considered as a complete convergence study, but only as 
a first test of the robustness of the results. Figure IB1I gives 
the stellar mass function, star formation density and cold gas 
mass function for the three boxes introduced in section [12] 
and for the same set of parameters given in table [4] 
The following conclusions can be drawn: 

(i) the model roughly converges in predicting the stellar 
content of bright galaxies, but does not for small galaxies. 
In particular, at increasing resolution (decreasing particle 
mass) the stellar mass function gets steeper at the low-mass 
end and, as a consequence, lowers at the knee. 

(ii) The star-formation density apparently converges at 
lower redshift, but gets increasingly large at z > 2. The main 
reason for the low-redshift convergence is that, to mimic the 
effect of reionization (following Benson et al. 2002; see sec- 
tion I2.6|l cooling is quenched in halos with circular veloc- 
ity larger than 50 km s _1 . The smallest sampled progenitor 
overtakes this limit at redshifts higher than 1, 2.8 and 7.5 
for the three runs, in order of decreasing particle mass. This 
means that each box is missing star-forming halos at red- 
shift higher than that limit, and explains why convergence 
in the star formation density is not visible in the figure. 

(iii) From the stellar mass function and star-formation 
rate, it is clear that the excess of small-mass galaxies no- 
ticeable for the highest resolution run is connected to the 
excess of star formation at high redshift, so that the exceed- 
ing dwarf galaxies will be very old objects. This unwanted 
feature is clearly connected to the excess of galaxies with 
stellar mass < 10 10 Mq at z ~ 1 found by Fontana et al. 
(2006) by comparing models (including MORGANA) to the 
results of the GOODS-MUSIC survey, and to the excess of 
small bulges visible in figure [7] Clearly, a source of feedback, 
which would limit the formation of dwarf galaxies at z > 1, 
is missing, and this lack is more and more evident when the 
particle mass decreases. 

(iv) With the low-V c cutoff, convergence is reached in 
the star- formation rate at z < 1, because only the run with 
the largest particle mass shows higher star formation density 
and a tail of galaxies more massive than 3 x 10 12 Mq . This is 
related to the efficiency of the quenching of the late cooling 
flows by radio- loud AGNs. This demonstrates that our self- 
consistent quenching needs a particle mass not larger than 
10 9 Mq. In case of poor mass resolution it is possible to 
apply a procedure similar to Croton et al. (2006) and Bower 
et al. (2006) to force the quenching of the cooling flow. In 
this case, whenever a cooling flow is present it is assumed 
that a fraction of that flow would immediately accrete onto 
the BH (subject to the Eddington rate) and release energy to 
the hot halo phase. Whenever this "fiducial" energy is higher 
than that lost by cooling, quenching is switched off. This 
"forced cooling" procedure is clearly less physical than that 
described in section [9721 and for a sufficiently small particle 
mass it gives a slightly stronger quenching of the cooling 
flows, but it can be considered as a rough numerical trick to 
achieve a good level of quenching when mass resolution is 
poor. This forced quenching procedure is used in Fontanot 
et al. (2006a). 

(v) The high-mass cutoff of the cold gas mass function is 




9 10 11 

Log M H , (M„) 



Figure Bl. Stellar mass function (upper panel; see figure 0, 
star formation density (mid panel; see figure[8]l and cold gas mass 
function (lower panel; see figure ITZt for galaxies in the three runs 
with particle mass 3.0 X 10 s M (blue dot-dashed lines), 1.0 X 10 9 
Mq (red dashed line) and 2.4 X 10 9 Mq (black continuous line). 
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the most sensitive prediction to mass resolution. We don't 
find any hint of convergence for this quantity at z = 0, 
despite the star content of the same galaxies is convergent. 
On the other hand, the low-mass tail of the same distribution 
is rather stable. 

We can then conclude that the model does not really 
converge with the resolution, and that convergence at large 
masses is obtained by using Benson et al.'s motivated recipe. 
Some relevant ingredient, able to hamper star formation in 
small-mass (but relatively high circular velocity) galaxies at 
high redshift, is still missing. 

This paper has been typeset from a Tp^X/ I^Tj^X file prepared 
by the author. 
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